Brain Synchronization and Multivariate Autoregressive (MVAR) Modeling in Cognitive Neurodynamics

This paper is a review of cognitive neurodynamics research as it pertains to recent advances in Multivariate Autoregressive (MVAR) modeling. Long-range synchronization between the frontoparietal network (FPN) and forebrain subcortical systems occurs when multiple neuronal actions are coordinated across time (Chafee and Goldman-Rakic, 1998), resulting in large-scale measurable activity in the EEG. This paper reviews the power and advantages of the MVAR method to analyze long-range synchronization between brain regions (Kaminski et al., 2016; Kaminski and Blinowska, 2017). It explores the synchronization expressed in neurocognitive networks that is observable in the local field potential (LFP), an EEG-like signal, and in fMRI time series. In recent years, the surge in MVAR modeling in cognitive neurodynamics experiments has highlighted the effectiveness of the method, particularly in analyzing continuous neural signals such as EEG and fMRI (Pereda et al., 2005). MVAR modeling has been particularly useful in identifying causality, a multichannel time-series measure that can only be consistently computed with multivariate processes. Due to the multivariate nature of neuronal communication, multiple non-linear multivariate-analysis models are successful, presenting results with much greater accuracy and speed than non-linear univariate-analysis methods. Granger’s framework provides causal information about neuronal flow using neural time and frequency analysis, comprising the basis of the MVAR model. Recent advancements in MVAR modeling have included Directed Transfer Function (DTF) and Partial Directed Coherence (PDC), multivariate methods based on MVAR modeling that are capable of determining causal influences and directed propagation of EEG activity. The related Granger causality is an increasingly popular tool for measuring directed functional interactions from time series data.

This paper is a review of cognitive neurodynamics research as it pertains to recent advances in Multivariate Autoregressive (MVAR) modeling. Longrange synchronization between the frontoparietal network (FPN) and forebrain subcortical systems occurs when multiple neuronal actions are coordinated across time (Chafee and Goldman-Rakic, 1998), resulting in large-scale measurable activity in the EEG. This paper reviews the power and advantages of the MVAR method to analyze long-range synchronization between brain regions (Kaminski et al., 2016;Kaminski and Blinowska, 2017). It explores the synchronization expressed in neurocognitive networks that is observable in the local field potential (LFP), an EEG-like signal, and in fMRI time series. In recent years, the surge in MVAR modeling in cognitive neurodynamics experiments has highlighted the effectiveness of the method, particularly in analyzing continuous neural signals such as EEG and fMRI (Pereda et al., 2005). MVAR modeling has been particularly useful in identifying causality, a multichannel time-series measure that can only be consistently computed with multivariate processes. Due to the multivariate nature of neuronal communication, multiple non-linear multivariate-analysis models are successful, presenting results with much greater accuracy and speed than non-linear univariate-analysis methods. Granger's framework provides causal information about neuronal flow using neural time and frequency analysis, comprising the basis of the MVAR model. Recent advancements in MVAR modeling have included Directed Transfer Function (DTF) and Partial Directed Coherence (PDC), multivariate methods based on MVAR modeling that are capable of determining causal influences and directed propagation of EEG activity. The related Granger causality is an increasingly popular tool for measuring directed functional interactions from time series data.

INTRODUCTION
Complex systems are of central interest to cognitive neuroscience. Any system composed of multiple moving parts is complex. The brain is a complex system. So is the cerebral (neo)cortex, that part of the brain most essential for cognition. Synchronization is a hallmark of complex systems. Brain synchronization within large-scale networks has been a primary research topic for many cognitive neuroscientists over the course of several decades. Large-scale coordination concerns brain synchronization within the broader context of coordinated structures in nature, a subject of study in physics, biophysics, and other disciplines. Directed functional connectivity is synchronization in regard to fMRI studies, and this paper discusses mathematical models for measuring directed functional connectivity. The mentioned directed functional connectivity is a key means of measuring synchronization within large-scale neocortical networks, or neurocognitive networks. Neurocognitive networks are defined as large-scale, distributed, interconnected systems of brain areas in the central nervous system that are joined together to perform a particular cognitive task (Bressler and Menon, 2010). Neurocognitive networks require a binding mechanism to function. The brain and cerebral cortex are complex systems that use binding mechanisms. They can be studied by applying directed functional connectivity methods because brain areas are linked by axon bundles, containing axons on which action potentials all travel in the same direction. This paper seeks to provide a review of Multivariate Autoregressive (MVAR) modeling in cognitive neurodynamics. Previous studies using MVAR modeling to analyze continuous neural signals are detailed, and further conclusions are drawn regarding given data and long-term implications.
Synchronization within neurocognitive networks concerns neural oscillations, brain activity that occurs in different brain areas that may or may not be spatially distant. Neural oscillations are unique in that they only display frequency, amplitude, and phase. In fact, when two or more neural oscillations are in synchronization, the exhibited frequency is always identical, a characteristic feature of the activity of neural oscillations. Amplitude and phase synchronization are both observed in the brain. Each is a possible mechanism of synchronization and may serve as a binding mechanism between neuronal populations within brain regions.
Considerable efforts have been made to model the directed functional connectivity in neurocognitive networks, and multiple mathematical models have been proposed to better interpret synchronization as a binding mechanism for neuronal populations in the brain. Of the models that have thus been developed, multivariate models are far more accurate and statistically significant than univariate non-linear models. Among other factors, the multi-dimensional nature of neuronal communication means that multivariate models are more ideal for interpreting directed functional connectivity and synchronization within the brain. Of the mathematical models that measure directed influences and causality within neurocognitive networks, Granger causality is the most basic and widely used, and is the most widely applicable to statistical contexts concerning causal influence within broader networks.
Granger's framework has been expanded on. Both partial directed coherence (PDC) and directed transfer function (DTF) exclusively interpret directed influences between time series in the multivariate framework (Faes et al., 2013). They both differ from Granger causality. PDC, a method introduced by Baccalá and Sameshima (2001), is an expansion of Granger causality that normalizes terms in the frequency domain by the total outflow at a site. DTF, a method introduced by Kaminski and Blinowska, normalizes frequency-domain terms by the total inflow at a site. It is an alternative causality model that normalizes directed influences by the sum of transfer functions entering the site, as opposed to the sum of transfer functions leaving a site that PDC measures. In most neural contexts, MVAR models using non-instantaneous effects are utilized to interpret rhythms. MVAR models, with and without instantaneous effects, are elaborated upon in the following sections.

Visuomotor Experiments in Macaque Monkeys
Macaque monkeys performed a visuomotor task with a GO/NO-GO response (Brovelli et al., 2004). Microelectrodes measured surface-to-depth Local Field Potentials (LFPs) from 15 sites across one hemisphere in the monkey's neocortex. During the task, the monkey depressed and held a lever in a random interval of 0.12 to 2.2 s. On GO trials, the monkey was rewarded with water if it released the lever within 500 msec after stimulus onset. Extracted data were preprocessed by selection and removal of trials containing contaminated or incorrect behavioral responses. Correct GOtrial LFP recordings were combined to result in approximately 900 trials per monkey. The data were then interpreted using MVAR spectral analysis (Geweke and Singleton, 1980;Geweke, 1982). Model coefficients were estimated by treating the analyzed LFP data as realizations of a common stochastic process.
Power, coherence, and relative phase spectral measures were estimated from the MVAR spectral matrix. The largest power peak by far was in the beta frequency range. Coherence and Granger causality spectra also contained prominent peaks in the beta frequency range. Coherence spectra were then analyzed for synchronized beta oscillations. A permutation distribution was created, and significance values were corrected using Dunn's multiple comparison procedure. Granger causality spectral analysis was also used to identify the predictability and relative strengths of influence at various cortical sites of the monkeys. Peak values and frequencies were identified and listed for each coherence and Granger causality spectrum at p < 0.005 (Figure 1).

FPN Synchronization During Working Memory in Macaque Monkey
Macaque monkeys performed a delayed match-to-sample task, while frontoparietal network (FPN) areas were monitored FIGURE 1 | Prestimulus beta-frequency coherence and (conditional spectral) Granger causality maps derived from the sensorimotor cortices of two monkeys (A,B). A self-generated hand press (-1250 to -500 msec) cued the monkey that a visual stimulus (0 msec onset; 100 msec duration) was soon to appear on a visual display screen. The stimulus was subsequently perceptually discriminated as part of a visual pattern discrimination task (water reward given on response trials 500 msec after stimulus onset) (Bressler et al., 1993). In each case, the pattern of synchronization (coherence) and Granger causality of beta-band oscillations from primary and secondary somatosensory and motor cortices consistent with execution of the hand press cue: somatosensory input is fed to the primary somatosensory cortex and motor output is transmitted from the primary motor cortex to the motor spinal cord to the hand muscles to execute the hand press. Sulci: Cs -Central; As -Anterior; Ls -Lunate; STs -Superior Temporal; IPs -IntraParietal (Figure modified from Brovelli et al., 2004). (Salazar et al., 2012). Two macaque monkeys matched the identity of a sample object, while recordings were made of broadband neuronal activity from 6 prefrontal cortical (PFC) and 6 posterior parietal cortical (PPC) sites (Figure 2). LFPs were recorded from monkeys A (over 27 days) and B (over 47 days). Time-frequency coherence and Granger causality spectra were computed for all FPN pairs (Figure 2). Pairs having significant spectra were identified, and peak values and frequencies were listed. Identity selectivity was also identified at each stimulus location.
Pairs with significant coherence selectivity index (CSI) were identified, and pairs having common CSI during the delay were grouped. Mean and variance of relative phase and power of betarange distributions were displayed. To determine the relation between cortical regions and FPN synchronization in working memory activity, fronto-parietal pairs showing significant CSIs were sorted. To determine cortical regions showing dominant activity in working memory, Granger causality evaluated the prediction that FPN synchronization is governed by synaptic influences in the PFC, originating in the PPC (Chauvette et al., 2012)

fMRI Blood-Oxygen-Level-Dependent Studies of Top-Down Influence in Human Neocortex
Human subjects performed a visual anticipation task designed to test the hypothesis that prefrontal and parietal cortices contain control areas that send top-down signals to visual cortex to modulate activity there and instantiate visual anticipatory attention (Bressler et al., 2011; Figure 3).

fMRI Blood-Oxygen-Level-Dependent Studies in the Human Dorsal Anterior Cingulate Cortex
Eleven adolescents participated in the study under proper ethical guidelines (Asemi et al., 2015; Figure 4). Participants were instructed to tap their right-hand forefinger as quickly as possible in response to a white flashing visual stimulus. fMRI Blood-Oxygen-Level-Dependent (BOLD) data were collected continuously across the study's conditions. A complete BOLD scan lasted 6.83 s. Functional volumetric images were preprocessed with SPM5 under standard protocol (Friston et al., 1995). Regions of Interest (ROIs) were dorsal anterior cingulate cortex (dACC), Supplementary Motor Area (SMA) and primary motor cortex (M1) in the left hemisphere. Eigenvariate time series were extracted from voxels centered in these ROIs from an effects of interest contrast. Values were averaged and the average eigenvalue time series underwent analysis. ROI time series were preprocessed by z-value normalization and outlier rejection.
Pearson product-moment correlation coefficients were computed between the dACC fMRI BOLD time series and that of SMA and M1 for each participant. DFC was then estimated from MVAR models. Correlations were compared between task and rest by the parametric paired t-test and the non-parametric Mann-Whitney-Wilcoxon signed-rank test.

Granger Causality
The general unrestricted equation of Granger causality is as follows: Granger causality is a popular multivariate model that provides causal interactions from time-series data. Granger causality is the main mathematical method by which dependencies on multiple variables can be investigated as opposed to undirected influences (Seth et al., 2015). Causal influences are estimated using the conditional Granger causality index (CGCI). Granger causality is often used in conjunction with Dynamic Causal Modeling (DCM), a non-linear system that utilizes the Bayesian Model Comparison. Despite early controversies regarding the full potential to capture causal mechanisms within functional neural circuits, Granger causality has been largely characterized as an accurate and reliable method to measure directed influence among the neural circuits that underly behavior, cognition, and perception. A recent extension of Granger causality is proving popular (West et al., 2020).

Multivariate Autoregressive Model Without Instantaneous Effects
The Multivariate Autoregressive Model without Instantaneous Effects is described by the following equation: X (n) = p k=1 A(k)X(n − k) + U(n), where p is the model order, A(k), k = 1,..., p, are M × M matrices contain the elements A ij (k) (Erla et al., 2009). MVAR models without instantaneous effects describe causality influences with lag and are generally less popular than MVAR models with instantaneous effects.

Multivariate Autoregressive Model With Instantaneous Effects
The Multivariate Autoregressive Model with Instantaneous Effects is described by the following equation: X (n) = FIGURE 3 | Top row. Top-down (top left) and bottom-up (top middle) Granger causality F-statistic histograms for a representative ROI pair, right aIPS and right V3A, in one subject. The critical value of F is 3.87 for significance (p < 0.05) in both directions. A larger fraction of the total number of voxel pairs (1064) had significant F statistics in the top-down (16.9%) than in the bottom-up (8.7%) direction. The schematic diagram (top right) showed Granger causality represented as arrows in the two directions between these ROIs on a standard brain image. Arrow thickness corresponds to the significant fraction, representing Granger causality strength. Bottom Top-down Granger causality before correct versus incorrect performance. (A) Grid specifying the significance of correct versus incorrect performance difference in a group analysis of top-down Granger causality: *p < 0.05, **p < 1.0 × 10 10 , ***p < 1.0 × 10 30 ; white, not significant. Each cell represents Granger causality from the row-labeled ROI to the column-labeled ROI. (B) The fraction of significant top-down Granger causality from right aIPS to left VP was significantly greater before correct (blue) than incorrect (red) performance in five of six subjects [Figure modified from Bressler (2018)].
Frontiers in Systems Neuroscience | www.frontiersin.org FIGURE 4 | Top The group mean and standard error of influence in both directions between the 2 ROIs (SMA to dACC, dACC to SMA) for Task (blue) and Rest (red) conditions. Post hoc paired t-tests revealed that the influence from dACC to SMA was significantly greater for Task than Rest, but that the influence from SMA to dACC did not significantly differ (**p < 0.01). Bottom The distribution across subjects for both directions between the 2 ROIs (SMA to dACC, dACC to SMA) and for Task (blue) and Rest (red) conditions, as box plots. Note that the distribution for the influence from dACC to SMA during the Task is both more compact than and elevated above the distribution during Rest (Figure modified from Asemi et al., 2015). p k=0 B k X n − k + W(n). The input noise W(n) = [W1(n),..., Wm(n)]ˆT is a vector of zero-mean uncorrelated processes with diagonal covariance matrix W (Erla et al., 2009). MVAR models with instantaneous effects describe causality effects with the instantaneous terms included, as opposed to MVAR models without instantaneous effects. MVAR models with instantaneous effects are generally more popular than MVAR models without instantaneous effects.

Partial Directed Coherence Derived From Granger Causality
Partial Directed Coherence is an MVAR model proposed by Baccalá and Sameshima (2001) that comprises a frequencydomain method for Granger causality. The Partial Directed Coherence from channel j to i is represented by the following equation: (Baccalá and Sameshima, 2001;Sameshima and Baccalá, 2014). Partial Directed Coherence displays directed influences between channels, normalized by the outflows from the j-th site.

Directed Transfer Function Derived From Granger Causality
Directed Transfer Function (DTF) is an MVAR model proposed by Kaminski and Blinowska (Kaminski and Blinowska, 1991;Kus et al., 2004) that is closely related to the PDC model. The DTF from channel j to i is represented by the following equation: (Kaminski and Blinowska, 1991). DTF is normalized by the sum of inflows to the i-th site. The model is applicable to a wide variety of neurophysiological contexts. To distinguish direct from indirect influences, direct Directed Transfer Function (dDTF) is utilized. Direct Directed Transfer Function from channel j to channel i is represented by the following equation: dDTF ji (f ) = F 2 ij C 2 ij (f ) (Korzeniewska et al., 2003). dDTF has been used, for example, to identify relations between primary motor sites.

Inherent Mechanisms of Multivariate Autoregressive Models
The MVAR models discussed are subsets of directed frequencydomain influences, under the general set of functional connectivity metrics. Functional connectivity is regarded as the incidence of statistically related neurophysiological events occurring in spatially distant regions of the brain (Pourahmadi and Noorbaloochi, 2016). Ongoing research on functional connectivity seeks to uncover large-scale neural circuits responsible for executive function and behavior. Modeling of brain activity is an increasingly critical undertaking in order to fully understand the mechanisms of brain function and connectivity. Although there have been considerable efforts to investigate and understand brain rhythms, future research will be needed to fully understand the mechanisms of brain activity, and in particular, calculated mean-power-spectra-displayed beta oscillatory activity at peaks surrounding 20 Hz, the oscillations being apparent in large-scale cortical networks. While classical univariate analysis is useful in particular contexts, multivariate time series analysis more fully captures the dynamic nature of brain connectivity and communication. The nature of neuronal connectivity is such that multiple inputs over a specified frequency and time-domain must be considered, thus justifying the need for MVAR modeling.
The methods of the discussed research are mainly mathematical. The following Multivariate Autoregressive Model is used to analyze continuous neural signals in the form of EEG and fMRI: X i,t = a i, 1, 1 X 1, t−1 + a i,1,2 X 1, t−2 + ... + a i,1,m X 1,t−m + a i,2,1 X 2,t−1 + a i,2,2 X 2,t−2 + ... + a i,2,m X 2,t−m + ... + a i,p,1 X p,t−1 + a i,p,2 X p,t−2 + ... + a i,p,m X p,t−m + e i,t . This equation can be further stated in matrix form as the following: where X t = [X 1t , X 2t , ..., X pt ] ∧ T are p data channels, m is the model order, A k are p × p coefficient matrices, and E t is the white noise residual error process vector. Spectral Analysis is then conducted from MVAR modeling using the Spectral Matrix

Visuomotor Experiments in Macaque Monkeys
The calculated BOLD mean power spectra displayed beta oscillatory activity at peaks surrounding 20 Hz mean. Coherence and Granger causality spectra also displayed peaks near 20 Hz. Coherence spectra displayed synchronization of oscillations.
Granger causal (GC) influences between cortical sites were mediated by beta range oscillations, as illustrated by the pronounced peaks surrounding 20 Hz in GC spectra. Coherence spectra were tested to identify significant beta peaks, a classic sign of synchronized beta LFP oscillations. Ultimately, coherence values were largest in three brain regions: the primary motor cortex, primary somatosensory cortex, and an area inferior to the intraparietal sulcus. GC oscillations within the beta oscillatory network were determined. GC spectra and coherence were closely linked, although GC influences were not found to be linked with the time delay values from phase spectra. Inferior posterior parietal sites exerted greater influence on the GC relations of primary motor sites.

Frontoparietal Network Synchronization During Visual Working Memory in Macaque Monkeys
The calculated CSI and WGC values demonstrated that FPN synchronization during visual working memory tasks is widespread, content-specific, and task-dependent during the delay period. PPC influences dominate the synchronization patterns (Verhoef et al., 2011). The findings of this experiment are generally consistent with past studies concerning the spatial attention modulation of inter-areal coherence (Verhoef et al., 2011). Short-term memories are represented as patterns that are specific to stimuli and are characteristic of synchronized activity widely distributed throughout the FPN (Tallon-Baudry et al., 2001). Ultimately, it was concluded that the FPN exerts top-down control of and transmits behavioral influences to the visual cortex during the delay period of the given task (Bressler and Richter, 2015;Richter et al., 2018). The beta frequency band was observed to contain the main spectral peak in the delay period, and the beta GC from the frontal cortex was especially strong during the delay period.
fMRI Blood-Oxygen-Level-Dependent Studies in the Human Neocortex Top-down influences from the frontal eye field and the intraparietal cortex to the visual cortex were found in a hemisphere of the human brain prior to an expected visual stimulus in relation to anticipatory visual spatial attention, suggesting that frontal and parietal control signals modulate sensory cortex (Bressler et al., 2011). Other site pairs showed bottom-up, but not top-down, influences.
fMRI Blood-Oxygen-Level-Dependent Studies in the Human Dorsal Anterior Cingulate Cortex Response data were not driven by frequency or periodicity, as covariance measures with frequency and periodicity did not have statistical significance. A statistically significant relationship was observed during Pseudo-random epochs for the relationship between the percentage of missed responses and age, suggesting that the frequency of missed responses decreases with age. Measured correlations between dACC and SMA, and between dACC and M1, were greater in the Periodic condition than the Pseudo-random condition by both conducted paired t-tests and signed-ranked tests. This finding suggests that dACC activity is tightly related to both areas when predictable stimuli are present. Ultimately, results demonstrated greater DFC in the top-down direction than the bottom-up direction, indicating the presence of top-down motor control, a major component of many complex behaviors, exerted by the dACC in this task. The employed mathematical models indicate with statistical significance that top-down synchronization exists in human brains. Furthermore, the results indicate that MVAR models are usable with statistical significance and accuracy to determine causal interactions among neuronal groups in neurocognitive networks. Brain synchronization has been found to be a binding mechanism within neurocognitive networks, as reaffirmed by past MVAR modeling studies in cognitive neuroscience research.

DISCUSSION
All four experiments tested the concept of brain synchronization as a binding mechanism within neurocognitive networks and provided evidence for top-down synchronization in monkey and human brains. Causal influences also suggested that top-down synchronization from the frontal cortex may compose the basis for predictive coding and selective attentional set in the brain. Isolation of sensory and motor areas in neurocognitive networks suggests that top-down FPN synchronization may contribute to human understanding and comprehension. The methods of MVAR modeling and GC are generally effective in measuring large-scale networks. Directional influences revealed by GC determine the electrical impulses of brain activity, especially when analyzing synchronous neural circuits and large-scale neuronal activity. These methods are particularly accurate when exogenous inputs direct intracerebral responses (Chang et al., 2012) but are widely applicable in all neurophysiological contexts. Modeling between cortical signals is necessary when determining electrical activity across cortical regions, thus explaining why univariate models are not as effective as multivariate methods are. MVAR modeling has provided evidence for synchronization and coordination in neurocognitive networks, as supported by numerous studies. Brain synchronization as a binding mechanism within neurocognitive networks continues to be heavily researched, and MVAR modeling has increased understanding of brain synchronization since the inception of these methods in neuroscience. MVAR models continue to evolve into an ever-expanding and more relevant method set of tools for analyzing neuronal communication, with new models being developed constantly. Despite recent advances in MVAR modeling (Pagnotta and Plomp, 2018), there remain numerous opportunities for growth in MVAR modeling of neuronal populations. First, Granger non-causality analysis depends on available model covariates, and thus is slightly limited in the method's scope when analyzing causality large-scale networks (Pourahmadi and Noorbaloochi, 2016). Second, the basic GC only models linear functions, once again restricting the scientific scope of this method (Seth et al., 2015). Since most normal EEGs are locally linear, this is not a great problem. Although these are minor hinderances in the applicability of GC, on which most MVAR models are based, there rests an even greater potential for analyzing synchronization and neural communication should these hinderances be removed.
The concept of brain synchronization as a binding mechanism within neurocognitive networks continues to be studied in great detail. Research continues to uncover new aspects of synchronization and coordination within the human brain.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by FAU HSC. The patients/participants provided their written informed consent to participate in this study. The animal study was reviewed and approved by FAU IACUC.

AUTHOR CONTRIBUTIONS
AK wrote the initial draft. SB conceptualized the review and edited the manuscript. All authors contributed to the article and approved the submitted version.