Introduction

The primate brain is highly interconnected between and within brain areas. This includes areas involved in sensory processing ¹. Strikingly, most computational models of brain activity in response to external natural stimuli do not take the recurrent architecture of brain networks into account. “Encoding” models often rely on simple input/output relationships such as general linear models in fMRI ², or temporal response functions in EEG/MEG ³. Interactions between brain areas are captured often just as instantaneous linear correlations that are referred to as “functional connectivity” when analyzing fMRI activity ⁴. Others capture synchronous activity in different brain areas by measuring phase locking of electrical neural signals ⁵. However, these measures of instantaneous correlation do not capture time delays inherent in recurrent dynamics. By taking temporal precedence into account with recurrent models the “Granger-causality” formalism can establish directed “connectivity”. This has been used to analyze both fMRI and electrical activity ^6–11.

The concept of functional connectivity was first developed to analyze neural activity during rest, where there are no obvious external signals to stimulate brain activity. But it is now also often used during passive exposure to a stimulus, such as watching movies ^12–15. A general observation of these studies is that a portion of the functional connectivity is preserved between rest and stimulus conditions, while some aspects are altered by the perceptual task, e.g. ^12,16. This should be no surprise, given that an external stimulus can drive multiple brain areas and thus induce correlations between these areas ¹⁷. Removing such stimulus-induced correlations by controlling for a common cause is standard practice in statistical modeling and causal inference ¹⁸. However, in studies that focus on functional connectivity in neuroscience, stimulus-induced correlations are often ignored when analyzing the correlation structure of neural signals. A notable exception is “dynamic causal modeling” ¹⁹. In this modeling approach the “input” can modulate functional connectivity. This is particularly important in the context of active behavioral tasks, where the common finding is that correlation structure changes with task states ²⁰.

In this study we are interested in “passive” tasks such as rest and movie watching. We will ask here whether, after removing stimulus-induced correlations, the intrinsic dynamic itself is preserved. Attempts to factor out the effects of the stimulus come from work on response variability. For instance, fMRI shows that variability across trials in motor cortex is due to an intrinsic “noise” which is linearly superimposed on a more reliable response to a simple motor action ²¹. Stimulus-response variability in the visual cortex has been attributed to variability of ongoing dynamic ^22,23. Some studies of electrical recordings from the visual cortex show that correlations of spiking activity between different recording locations are largely unaffected by visual stimulation ²⁴. Yet, other studies show that visual input affects local correlation in the visual cortex ^25–27 and across the brain ²⁸.

The technical challenge when addressing these questions is to separate the direct effect of the stimulus from the intrinsic recurrent dynamic. Here we propose to separate these effects by modeling them simultaneously with the simplest possible model, namely, linear intrinsic effects between brain areas and linear responses to extrinsic input. A mathematical model that implements this is the vector-autoregressive model with external input (VARX). This model is well established in the field of linear systems ²⁹ and econometrics ³⁰, where it is used to capture intrinsic dynamics in the presence of an external input. The VARX model is an extension of the VAR model that is routinely used to establish “Granger-causality” in neuroscience (cited above). In the VARX model, Granger analysis provides a measure of statistical significance for the intrinsic dynamic as well as the external input, in addition to directionality for the intrinsic effects, all as part of a single model ³¹.

While linear systems are an inadequate model of neuronal dynamics, they remain an important tool to understand neural representations because of their conceptual simplicity. They are routinely used for event-related fMRI analysis but also for “encoding models” to link non-linear features of continuous stimuli to neural responses. They have been used to analyze responses to video in fMRI ³², to speech in EEG ³³ or to audio in intracranial EEG ³⁴. They are even used to analyze the encoding in deep-neural network models ³⁵. Here we use a classic linear model to combine two canonical concepts in neuroscience, which have thus far remained separated, namely, that of “encoding models” ³² and “functional connectivity” models ⁶. We will use this to analyze whole-brain, intracranial EEG in human subjects at rest, and while they watch videos. Our main finding is that the recurrent dynamic observed during rest is only minimally altered by watching videos. Instead, the brain’s response to naturalistic stimulus appears to be substantially shaped by the same endogenous dynamic of the brain observed during rest.

Methods

The vector-autoregressive model with external input (VARX) falls within a group of well-established linear models used in neuroscience (see Table 1). Prominent examples in this group are the generalized linear model (GLM), dynamic causal model (DCM) and temporal response functions (TRF). While these models have been extensively used for neural signal analysis, the VARX model has not. We start therefore with a brief introduction. For more details please refer to ³¹

Table 1:

VARX model

The VARX model explains a time-varying vectorial signal y(t) as the result of an intrinsic autoregressive feedback driven by an innovation process e(t) and an extrinsic (We adopt here the terminology of “Intrinsic” and “extrinsic” as it is commonly used in neuroscience and psychology. In system modeling and econometrics, where the VARX model is prevalent, the more common terminology is “endogenous” and “exogenous”, meaning effectively the same thing) input x (t):

A* and B* represent convolutions with appropriately sized matrices of causal filters with lengths n_a and n_b respectively. The innovation is assumed to be uncorrelated in time and has therefore a uniform spectrum. The recurrence in A modifies this spectrum to match the spectrum of y(t), thereby capturing the intrinsic recurrent dynamic. The filter B injects a filtered version of the extrinsic input x (t) into this recurrent dynamic. The role of each of these terms for brain activity is explained in Fig. 1.

Figure 1:

Filter matrices A and B are unknown and can be estimated from the observed history of x (t) and y(t)using ordinary least squares (OLS). The objective for the optimal model is to minimize the power of the unobserved innovation process e(t):

Granger analysis

The innovation is also the prediction error, for predicting y(t) from the past y(t−1) and input x (t). In the Granger formalism the prediction error is calculated with all predictors included (error of the full model, σ_f) or with individual dimension in y(t−1)or x (t)omitted (error of the reduced models, σ_r) ³⁶. To quantify the “effect” of the specific dimension one can take the ratio of these errors ³⁷ leading to the test statistic D known as the “deviance”. When the number of samples is T large, the deviance follows the Chi-square distribution with cumulative density F, from which one can compute a p-value:

The p-value quantifies the probability that a specific connection in A or B is zero. The “generalized” R ^{2 38} serves as a measure of effect size, capturing the strength of each connection (D, p and R can be computed for each connection in matrix A or B). While this Granger formalism is well established in the context of estimating, A i.e. VAR models, to our knowledge, it has not been used in the context of estimating B, i.e. VARX or TRF models.

Overall system response

The overall brain response to the stimulus for the VARX model is given by the system impulse response (written here in the z-domain, or Fourier domain):

What we see here is that the system response H is factorized into an autoregressive (AR) filter and A a moving average (MA) filter B. When modeled as a single MA filter, the total system response has been called the “multivariate Temporal Response Function” (mTRF) in the neuroscience community ³⁹. We found that the VARX estimate of H is nearly identical to the estimated mTRF ³¹. In other words, B and A are a valid factorization of the mTRF into immediate extrinsic versus recurrent intrinsic effects.

Note that the extrinsic effects captured with filters B are specific (every stimulus dimension has a specific effect on each brain area), whereas the endogenous dynamic propagates this initial effect to all connected brain areas via matrix, A effectively mixing and adding the responses of all stimulus dimensions. Therefore, this factorization separates stimulus-specific effects from the shared endogenous dynamic.

Relation to common neural signal models

The VARX model fits naturally into the existing family of models used for neural signals analysis (Table 1). While they differ in the formulation and statistical assumptions, their defining equations have a similar general form with the following attributes:

An important simplifying assumption for the mTRF, VAR, and VARX models is that y(t) is observable with additive normal distributed innovation. As a result, parameter estimation can use ordinary least squares, which is fast to compute. In contrast, GLM, DCM, and some variants of VAR models assume that y(t) is not directly observable, and needs to be estimated in addition to the unknown parameters A or B. The same is true for the basic “output error” model in linear systems theory ²⁹. This requires slower iterative algorithms, such as expectation maximization. As a result, these models are often limited to small networks (The original DCM proposed for fMRI included an added complication of modeling the hemodynamic response, which amounts to adding a temporal filter to each output node and prior to adding observation noise.) of a few nodes to test specific alternative hypotheses ⁴². In contrast, here we will analyze 100-200 channels per subject to draw general conclusions about overall brain organization.

Validation of connectivity estimate on whole-brain neural mass model

To validate the interpretation that A is a model of “connectivity”, we simulated neural activity for a whole-brain neural mass model ⁴³. We used the default model of the neurolib python library (“ALNModel”), which is a mean-field approximation of adaptive exponential integrate-and-fire neurons. This model can generate simulated mean firing rates in 80 brain areas based on connectivity and delay matrices determined with diffusion tensor imaging (DTI). We used 5 min of “resting state” activity (no added stimulus, simulated at 0.1ms resolution, subsequently downsampled to 100Hz). The true connectivity matrix from DTI (Fig. 2A) appears to be similar to the effect size estimate R for the recurrent connections A in the VARX model with no input (Fig. 2B). Following ⁴⁴ we compare the two as a scatter plot (Fig. 2C) and observed a Spearman correlation of 0.69. For comparison, we also used the sparse-inverse covariance method to recover structural connectivity from functional connectivity. This method is more sensitive than others in detecting network connections ⁴⁵ and uses the graphical lasso algorithm ⁴⁶. The resulting connectivity estimate (Fig. 2D) only achieves a Spearman correlation of 0.52. We note that the structural connectivity determined with DTI is largely symmetric. When enhancing the asymmetry the VARX model is not as accurate, but correctly recovers the direction of the asymmetry (Fig. S1).

Figure 2:

Intracranial EEG recordings and stimulus features

We analyzed intracranial EEG and simultaneous eye-tracking data recorded from patients (N=21, mean age 37.81 years, age range 19-58 years, 9 female, Table S1) during rest and while they watched various video clips. Three patients underwent two implantations and recordings at different times resulting in a total of 24 recording sessions with a total of 4,962 recording channels. The video clips included animations with speech (‘Despicable Me’, two different clips, 10 min each, in English and Hungarian), an animated short film with a mostly visual narrative and music, shown twice (‘The Present’, 4.3 min), and three clips of documentaries of macaques (‘Monkey’, 5 min each, without sound) ⁴⁸. In addition to the clips from the previous analysis, we included a movie clip of abstract animations (‘Inscapes’, 10 min) ⁴⁹, and an eyes-open resting state with maintained fixation (‘Resting state’, 5 min). In total, we recorded up to 59.3 minutes of data for each patient (Table S1). Two patients did not complete both movie watching and resting state (Pat_5 & Pat_16) and were not included in the analysis that compares the two conditions.

Neural signals were preprocessed as previously described to reduce noise ⁴⁸. We re-reference signals in a bipolar montage to ensure analysis of local activity. We analyze local field potentials (LFPs) and broadband high-frequency activity (BHA) power. BHA is the power of the signal bandpass filtered between 70-150Hz. We perform analysis on both signals after downsampling to 60Hz. Example traces of y(t) for LFP and BHA are shown in Fig. 1B&C.

We extract three features of the movies that serve as external inputs for the VARX model: fixation onset, film cuts and sound envelope (Fig. 3G). Fixation onset and film cuts are represented in x(t) as pulse trains with pulses occurring at the time of these events ⁴⁸. Sound envelope is computed as the absolute value of the Hilbert transform of the sound from the movie files and varies continuously. The envelope is downsampled to 60 Hz. All videos and resting state include fixations. The video ‘Inscapes’ and resting state do not include film cuts as external input. The ‘Monkey’ video clips and resting state do not include the sound envelope as input features, but do include fixation onsets. When a feature is not available it is replaced with features from a different recording. Therefore, the statistics of the feature are consistent, but not aligned to the neural recording. When comparing models with different features we always keep the number of input variables consistent between models to avoid a bias by the number of free parameters of the model. Features that are not considered in the analysis are shuffled in time by a circular shift by half the duration of the signals.

Figure 3:

The VARX models were fitted to data with the matlab version of the code ³¹. For all analyses we use filters of 600 ms length for inputs (n_b=36 samples for VARX models, L=36 samples for mTRFmodels). Delays for connections between channels are set to 100ms (n_a=6 samples) for both LFP and BHA signals. Increasing the number of delays n_a, increases estimated effect size R, however, larger values lead to overfitting, i.e. less significant connections (Fig. S1). Values around n_a=6 samples achieve a balance between goodness of fit and overfitting (Fig. S2). The regularization parameter was set to λ=0.3.

Connectivity plots are created with nilearn’s plot_connectome() function (Fig. 4) ⁵⁰. We plot only significant connections (p<0.001). Surface plots of T1w/T2w ratios and directionality of connections are created using the field-echos repository ^51,52. T1wT2w maps ⁵³ are obtained from the neuromaps repository ^54,55, and transformed to the freesurfer surface using the fslr_to_fsaverage() function ^56,57.

Figure 4:

The length of responses for Fig. 5 is computed as the ‘peak widths’ argument of Matlab’s findpeaks() function. Power is computed as the average of the instantaneous power, i.e. the square of the weight at each delay of the filters.

Figure 5:

Data and code availability

The raw data reported in this study cannot be deposited in a public repository because of patient privacy concerns. To request access, contact The Feinstein Institutes for Medical Research, through Dr. Stephan Bickel. In addition, processed datasets derived from these data have been deposited at https://doi.org/10.17605/OSF.IO/VC25T and are publicly available as of the date of publication.

All original code has been deposited at https://github.com/MaxNentwich/varx_demo.

Results

Extrinsic input leads to spurious intrinsic connectivity

To determine the effect of the extrinsic inputs on connectivity estimates we either fit a VARX model or a VAR model (i.e. a VARX model with no external input). We analyze LFP data on all available recordings, movies and resting state for all N=24 recording sessions. As extrinsic inputs we included film cuts, fixation onset, and sound envelope. VAR models contain the same external inputs as the VARX model, but the time alignment is disrupted by a circular shuffle. This keeps the number of parameters in different models constant and ensures the inputs have the same covariance structure. We found a similar connectivity structure for the estimated VAR and VARX models (Fig. 3A and 3B). However, they vary systematically in the number of significant recurrent connections A (those with p<.0001, Fig. 3D), which drops when adding inputs (median=-8.7*10⁻⁴, p<.0001, N=24, Wicoxon). The effect sizes R also significantly decreases in the VARX model (Fig. 3E, median=-1.9*10⁻⁵, p<.0001, N=24, Wicoxon). Therefore, accounting for the external input removes spurious “connections”. We also analyzed how much each of these inputs contributed to this effect (Fig. 3F). Out of the three input features considered, models including fixations and cuts decrease effect size more than models with sound envelope (fixations vs. sound, medianΔR=-1.0*10⁻⁵, p<.0001, N=24; cuts vs. sound: median ΔR=-3.8*10⁻⁶, p<.0001, N=24; Wilcoxon, uncorrected). The model including the combination of all three features has a smaller effect size R for A than models with any individual input feature (all vs. fixations: median ΔR=-6.5*10⁻⁶, p<.0001, N=24; all vs. cuts: medianΔR=-1.2*10⁻⁵, p<.0001, N=24; all vs. sound: medianΔR=-1.9*10⁻⁵, p<.0001, N=24; Wilcoxon, uncorrected). Thus, adding more input features further reduces the strength of intrinsic “connections”. These results are also reflected in the analysis of BHA signals (Fig. S3).

Recurrent connectivity unchanged during movies and rest

Next we compared intrinsic “connectivity” between movie watching and rest (Fig. 4A-D). In the rest condition subjects have a fixation cross on a gray background. This obviously reduces the size and number of saccades as compared to movie watching, but does not abolish them (Fig. S4). We therefore use a VARX model including fixation onset as extrinsic variable in both cases. Movies include film cuts and the sound envelope as external inputs. To control for the number of free parameters, we include copies of the film cut and sound envelope features from the movies to the resting state model. Remarkably, the number of significant recurrent connections in A were not detectably different between movie watching and rest (Fig. 4C, median=-0.0019, p=0.19, N=22, Wilcoxon), as is the effect size (Fig. 4D, median=-9*10⁻⁵, p=0.14, N=22, Wilcoxon). One caveat to this conclusion is that the signal we analyzed was only 5 minutes long for the movie and rest conditions, and longer records may have revealed small differences. However, even on 5 minutes of data we observe a decrease in R values when including external inputs (Fig. S5). Connectivity of BHA between movie and rest does also not differ significantly (Fig. S6). Using different segments of movies, in some cases we find a small reduction of significant connections in movie watching compared to resting state conditions (Fig. S7). However, overall, differences in the intrinsic connectivity between movie and rest, if they exist, are less systematic than the effect of the stimulus.

Recurrent dynamic enhances and prolongs stimulus responses

We also compared the immediate exogenous effect B with the total system response H, which includes the additional effect of the recurrent dynamic A. We estimate B with the VARX model (Fig. 5A) on data during video watching, and estimate the total response H directly using temporal response functions (Fig. 5B). Both models include fixation onset, film cuts and sound envelope as external inputs. We compare the power and length of filters from both models (Fig. 5C-D). We compare responses in channels with significant effects of B (FDR correction, α=0.05). We see that the total response fixation onset is significantly stronger (Fig. 5C, medianΔ=-5.4*10⁻⁵, p<.0001, N=23, Wilcoxon) and longer than the immediate effect B (Fig. 5D, medianΔ=-21.72ms, p<.0001, N=23, Wilcoxon). The same effect is observed for other input features and for BHA responses (Fig. S8). This suggests that the total response of the brain to these external inputs is dominated by the recurrent dynamic of the brain.

Results are similar for VARX models of BHA and LFP

We repeated the same analyses of Figures 3-5 with broadband high frequency activity (BHA). While LFP are thought to capture dendritic currents, BHA is correlated with neuronal firing rates in the vicinity of an electrode. Generally we find a more sparse recurrent connectivity for BHA as compared to LFP (compare Fig. 3&4 with Fig. S3&S6). Perhaps this is expected, given that LFP covers a broader frequency range. Regardless of this overall difference, we find similar results when analyzing BHA with the VARX model. Namely, taking the extrinsic input into account removed stimulus-induced intrinsic “connections” (Fig. S3); the resulting model of the recurrent dynamic is indistinguishable between watching movies and rest (Fig. S6); and responses to the stimulus are stronger and more prolonged when separately modeling the effect of recurrent connectivity (Fig. S8). In the Discussion section we will argue that some of these results are expected in general when decomposing the total system response into extrinsic and intrinsic effects. What we did not necessarily expect is that the intrinsic dynamics is similar during movies and rest for both LFP and BHA.

Intrinsic “noise” in BHA is reduced by external stimulus

So far we have discussed the mean response captured by B and the recurrent activity mediated by A. We now want to analyze whether the external input modulates the variability of the internal dynamic. As a metric of internal variability we measured the power of the intrinsic innovation process e(t). For the LFP signal we see a drop in power during movies as compared to rest, for both the original signal y(t) (Fig. S10A) and the model’s innovation process e(t) (Fig. S10B). Notable is the stronger oscillatory activity during rest (Fig. S10A). In this example we see a drop in power in the theta/alpha band (5-11 Hz) during movie watching across all electrodes (Fig. S10A, dotter lines). We observe similar narrow-band drop in power in most patients, albeit at different frequencies (not shown). When analyzing BHA, we find no difference in power of the innovation process between movie and rest, but we do find a drop in power relative to the overall BHA signals for some channels (Fig. 6B). These channels seem to coincide with channels that responded to the external stimuli, i.e. channels with a significant effect in B (Fig. 6A). If we take for each subject the median relative power for responsive channels (median among those with p<0.0001), then we find that relative power drops for nearly all subjects (Fig 6D, Wilcoxon rank sum test, p=2.6e-06, N=21). The motivation for analyzing only responsive channels comes from a simple gain adaptation (Fig. S11). Gain adaptation keeps the power of y(t) constant, so that the extra power injected by the stimulus implicitly reduces the relative power of the innovation process. This effect is specific to channels receiving external input (Figl. S11D) and absent in a linear system without gain adaptation (Fig. S11C). To demonstrate that this simple gain adaptation can explain the noise quenching in the neural data, we simulated data with the gain adaptation model (Fig. 6C) using parameters estimated for the example subject of Fig. 6A/B.

Figure 6:

Direction of connectivity differs with cortical hierarchy

Finally, we measured the directionality of the recurrent connections in the LFPs by analyzing the structure of the resulting matrices R of all subjects. Columns in R represent outgoing connections, while rows are incoming connections. Therefore, the difference of R − R ^T (Fig. 7A) averaged along a column has positive values if a node has overall stronger outgoing connections, and negative values if it has stronger incoming connections. We measured this directionality for each channel across all subjects and averaged also across channels within parcels of the Desikan-Killiany atlas (N=35 regions of interest, Fig. 7B) ⁵⁸. We expected this to co-vary with “cortical hierarchy”. To test this, we compared this asymmetry metric with the T1w/T2w ratio, which captures gray matter myelination and is used as an indirect measure of cortical hierarchy ^51,59. We also average T1w/T2w ratio in the same parcels of the Desikan-Killiany atlas (Fig. 7B). Cortical areas showing more outgoing connections (R − R ^T > 0) have lower T1w/T2w ratio, which are located higher on the cortical hierarchy (Pearson’s r = 0.39, p = 0.023, Fig. 7C). BHA analysis shows the same trend (Fig. S12).

Figure 7:

Discussion

Our results suggest that intrinsic dynamics are not substantially altered during watching movies as compared to rest. Instead, the external stimulus reverberates in the recurrent network with the same dynamic as during rest. The duration and magnitude of response is in large part a result of this recurrent dynamic.

Response to extrinsic input versus intrinsic dynamics

Previous literature does often not distinguish between intrinsic connectivity and extrinsic effects. As a result, similarities and differences between rest and stimulus conditions reported previously, do not draw a firm conclusion as to whether “functional connectivity” is preserved, e.g. ^12,16. By systematically factoring out the effect of the external input we conclude here that the intrinsic dynamic is unaltered. If one fails to factor out the effect of the stimulus, one may mistake the stimulus-induced correlations for changes in “functional connectivity”.

In this work we focused on “passive” tasks, i.e. resting with gaze on a fixation point, versus watching movies without any associated tasks. We did not analyze data during an active task requiring behavioral responses. The literature on active tasks emphasizes “state change” in functional connectivity. ^14,20,60 Efforts to factor out task-evoked activity when computing functional connectivity concord with our conclusions that connectivity is inflated by a task ¹⁷. Nevertheless, we hesitate extrapolating our findings to active tasks, as we have not analyzed such data.

Conventional “encoding” models, such as temporal response functions, capture the total response H of the brain to an external stimulus. Here we factored this into a moving average filter B, followed by and autoregressive filter A. The important observation is that this intrinsic dynamic governed by A does not change during stimulus processing. Arguably then, the role of the initial responses B is to shape the input to be processed by the existing intrinsic dynamic. This interpretation is consistent with the view of “the brain from the inside out” advocated by György Buzsáki ⁶¹. In this view, learning of a stimulus representation consists in learning a mapping of the external stimulus to an existing intrinsic dynamic of the brain.

Similar findings for LFP and BHA

We found a more sparse recurrent connectivity for BHA as compared to LFP. This may be expected because correlations in lower frequencies (that dominate LFPs) reaches over longer distances compared to correlations in higher frequencies (e.g. Muller et al., 2016). BHA has been linked to a mixture of neuronal firing and dendritic currents ⁶², in contrast to LFP, which is thought to originate from widespread dendritic currents. Despite the observed differences in sparsity, for both LFP and BHA we found that modeling the recurrent dynamic removed spurious intrinsic connections. Removal of spurious effects when controlling for a common cause is a generic finding in multivariate statistical models. We also found for both LFP and BHA that the duration and strength of stimulus responses can be largely attributed to the recurrent dynamic. Arguably, this is a generic feature of an autoregressive model, as it more readily captures longer impulse responses. However, the extrinsic filters B in principle have an advantage as they can be fit to each stimulus and brain location. In contrast, the recurrent filters A are constrained by having to capture a shared dynamic for all stimulus dimensions. Thus, the predominance of the recurrent dynamic in the total system response is not a trivial result of the factorization into intrinsic and extrinsic effects. Finally, we did not necessarily expect that the intrinsic connectivity is preserved between movie and rest in both LFP and BHA. This consistency may be due to a variety of processes that are constant across conditions, such as internal thought, body and eye movements. Active sensing through eye movements, for example, influences activity in a global network ^63,64, and likely accounts for part of the common source of correlations across conditions.

Stimulus-induced reduction of noise in the intrinsic activity

One difference we did find between LFP and BHA is the intrinsic innovation process, i.e. the internal sources of variability or “noise”. For both BHA and LFP we saw a drop in the magnitude of signal fluctuations during the movie watching condition. For the BHA but not the LFP, this was explained as a drop in intrinsic noise. Specifically, for BHA there was less relative power in the intrinsic “noise” for channels that are responsive to the stimulus. This is consistent with the notion that response variability is due to variability of intrinsic activity ²² which is found to decrease across the brain with the onset of an external stimulus ⁶⁵. This type of noise quenching has been associated with increased attention ⁶⁶ and improved visual discrimination performance ⁶⁷. The effect we found here can be explained by a VARX model with the addition of a divisive gain adaptation mechanism that keeps the total power of brain activity constant. When the input injects additional power, this nonlinear gain adaptation implicitly reduces the contribution of the intrinsic noise to the total power.

We also observed an overall drop in LFP power during movie watching. This phenomenon was strongest in oscillatory bands, with frequencies in theta to beta band differing across subjects. In scalp EEG, noise quenching is associated with a similar overall drop in power with the stimulus ⁶⁶. This quenching of neural variability was also found to reduce correlation between brian areas for fMRI and neural spiking ²⁶. Both fMRI and neural spiking correlated with BHA ⁶⁸. This is at odds with our finding that intrinsic connectivity in BHA does not change significantly between movie and rest. However, we can not rule out such differences on longer recordings.

Stimulus features

During the movie and rest periods, we utilized fixation onset to capture activity that is time-locked to visual processing because subjects move their eyes even during rest. We also incorporated the sound envelope, a prominent feature known for capturing the dominant audio-induced variance in scalp EEG ³³. In addition, we included film cuts as features, as we had previously demonstrated that they dominate the response in the BHA across the brain ⁴⁸. While other basic visual features such as overall optic flow or fixations on faces elicited responses in the BHA, their contribution was relatively smaller. The analysis is not limited to these few features, and future research should explore which stimulus features capture variance in the data and how they affect the apparent intrinsic connectivity. There is a substantial body of literature on encoding models of semantic features, where nonlinear features of a continuous natural stimulus are extracted and then linearly regressed against fMRI ^69,70 or EEG ⁷¹. This work can be directly replicated with the VARX model which further models the intrinsic connectivity.

Alternative approaches

The traditional VAR model has been used extensively in neuroscience to establish directed “Granger causal” connections ⁴¹. This approach has been very fruitful and found numerous extensions, e.g. ^10,11. However, these model implementations do not specifically account for an external input.

A few methods have attempted to model the effect of varying task conditions on functional connectivity, mostly in the analysis of fMRI. One approach is to first model the task-evoked responses, equivalent to estimating B alone, and then compute the conventional “functional connectivity”, i.e. the correlation matrix, on the residuals e(t) ⁷². Others suggested to estimate B in multiple time windows and then estimate a “task related functional connectivity” by correlating the multiple B over time windows ⁷³. It is not clear that these ad-hoc methods systematically separate intrinsic from extrinsic factors.

A more principled modeling approach is “dynamic causal modeling” (DCM) ¹⁹ and extensions thereof ⁷⁴. Similar to the VARX model, DCM includes intrinsic and extrinsic effects A and B. However, the modeling is limited to first-order dynamics (i.e. n_a=n_a=1). Instead, the DCM includes a multiplicative interaction of extrinsic input x(t) on the connectivity A, which does not exist in the VARX model. This interaction has been used to explicitly model a change in intrinsic connectivity with task conditions. Here we found that this may not be necessary for intracranial EEG. A practical advantage of the VARX model is the assumption that the neural activity is directly observed. Instead, many existing models assume an error in the observations, which triggers computationally intensive estimation algorithms, typically the expectation maximization algorithm. The same is true for the “output error” model in linear systems theory ²⁹. As a result, these models are often limited to small networks (The original DCM proposed for fMRI included an added complication of modeling the hemodynamic response, which amounts to adding a temporal filter to each output node and prior to adding observation noise.) to test specific alternative hypotheses ⁴². In contrast, here we have analyzed 100-200 channels per subject across the brain, and have drawn more general conclusions about whole-brain activity.

Caveats

The lack of a significant difference in recurrent connectivity between stimulus and rest should be interpreted with care. As usual, lack of evidence is not evidence for the lack of an effect. We saw no change in the number of recurrent connections between movie and rest, either for the LFPs or BHA activity. However, in individual movie segments small differences were observed (Fig. S7). It is possible that regressing out a richer stimulus characterization would have removed additional stimulus-induced correlation, only enhancing this small difference between movie and rest. We were also limited to 5 minutes of data in the direct comparison of movie and resting state data. Longer recordings might further enhance differences. Higher recurrent connectivity in the LFP during rest would be consistent with the more synchronized state we saw in rest, as reflected by larger oscillatory activity.

We find a correlation of DTI structural connectivity used in a model with a VARX estimate of 0.70. That is considered a relatively large value compared to other studies that attempt to recover DTI connectivity from the correlation structure of fMRI activity ⁴⁴. A Caveat is that this was done on a biophysical model of firing rate, not fMRI, and we have not explored the parameters of the model that might affect the results.

We used fixation onsets as external input, but it should be noted that they are tightly correlated in time with saccade onsets (there is only about a 30 ms jitter between the two, depending on saccade amplitude). While saccades are driven by visual movement, they are generated by the brain itself and arguably could also be seen as intrinsic. The same is true for all motor behaviors, most of which cause a corresponding sensory response, similar to the visual response following a saccade. Including them as external input is a modeling choice we have made here, but it is important to acknowledge that fixation onsets can therefore have “acausal” components ⁴⁸. By “acausal” we mean a fixation-locked response that precedes the fixation onset and is due to the neural activity leading up to the saccade and subsequent fication. Such acausal responses can be captured by the VARX Granger formalism by delaying the input relative to the neural activity, which we have not done here.

The correlation between the average incoming and outgoing connections and cortical hierarchy (Fig. 7) is not significant when normalizing for the number of electrodes in each region of interest. Regions in the temporal lobe with a large number of electrodes might drive this correlation. A more fine grained analysis in these regions could be the goal of future analysis.

Conclusion

We analyzed whole-brain intracranial recordings in human subjects at rest and while they watched videos. We used a model that separates intrinsic dynamics from extrinsic effects. We found that the recurrent dynamic observed during rest is largely unaltered when watching movies. Instead, the brain’s response to the audiovisual stimuli appears to be substantially shaped by its endogenous dynamic. The reduction in intrinsic variance observed during an extrinsic stimulus may be the result of neuronal gain adaptation.

Acknowledgements

We would like to thank Chris Honey for advice on the model validation with simulations and related references. We like to thank Behtash Babadi for help on the development of the Granger formalism for the VARX model. This work was supported in part by the NIH through grants P50 MH109429 and R01DC019979.

Additional information

Author contributions

Conceptualization, L.C.P. and M.N.; Methodology, L.C.P.; Software, L.C.P. and M.N.; Formal Analysis: M.N and L.C.P.; Investigation, M.N. and M.L.; Resources, S.B.; Writing – Original Draft, L.C.P. and M.N.; Writing – Review & Editing, L.C.P., M.N., M.L., C.E.S., S.B.; Funding Acquisition, L.C.P., C.E.S. and S.B.; Visualization, M.N.; Supervision, L.C.P., C.E.S., S.B.

Additional files

Supplemental information