Multivariate stochastic volatility modeling of neural data
Abstract
Because multivariate autoregressive models have failed to adequately account for the complexity of neural signals, researchers have predominantly relied on nonparametric methods when studying the relations between brain and behavior. Using medial temporal lobe (MTL) recordings from 96 neurosurgical patients, we show that time series models with volatility described by a multivariate stochastic latentvariable process and lagged interactions between signals in different brain regions provide new insights into the dynamics of brain function. The implied volatility inferred from our process positively correlates with highfrequency spectral activity, a signal that correlates with neuronal activity. We show that volatility features derived from our model can reliably decode memory states, and that this classifier performs as well as those using spectral features. Using the directional connections between brain regions during complex cognitive process provided by the model, we uncovered perirhinalhippocampal desynchronization in the MTL regions that is associated with successful memory encoding.
https://doi.org/10.7554/eLife.42950.001Introduction
Recent advances in neuroscience have enabled researchers to measure brain function with both high spatial and temporal resolution, leading to significant advances in our ability to relate complex behaviors to underlying neural signals. Because neural activity gives rise to electrical potentials, much of our knowledge concerning the neural correlates of cognition derive from the analyses of multielectrode recordings, which yield a multivariate time series of voltage recorded at varying brain locations (denoted here as ${\text{\mathbf{y}}}_{t}$). Such signals may be measured noninvasively, using scalp electroencephalography (EEG), or invasively, using subdural grids or intraparenchymal depth electrodes in human neurosurgical patients. In recent years, intracranially recorded (iEEG) signals have yielded detailed information on correlations between timeseries measures and a wide range of behaviors including perception, attention, learning, memory, language, problem solving and decision making (Jacobs and Kahana, 2010).
Whereas other fields that grapple with complex multivariate time series have made effective use of parametric models such as economics and engineering (Kim et al., 1998; Blanchard and Simon, 2001; West, 1996), neuroscientists largely ceded early parametric approaches (e.g. linear autoregressive models) in favor of nonparametric spectral decomposition methods, as a means of uncovering features of neural activity that may correlate with behavior. A strength of these nonparametric methods is that they have enabled researchers to link fluctuations in iEEG signals to lowfrequency neural oscillations observed during certain behavioral or cognitive states, such as slowwave sleep (Landolt et al., 1996; Chauvette et al., 2011; Nir et al., 2011), eye closure (Klimesch, 1999; Goldman et al., 2002; Laufs et al., 2003; Barry et al., 2007) or spatial exploration (Kahana et al., 2001; Raghavachari et al., 2001; Caplan et al., 2003; Ekstrom et al., 2005; Byrne et al., 2007). Highfrequency neural activity, which has also been linked to a variety of cognitive and behavioral states (Maloney et al., 1997; Herrmann et al., 2004; Canolty et al., 2006), is less clearly oscillatory, and may reflect asynchronous stochastic volatility of the underlying EEG signal (Burke et al., 2015).
Although spectral analysis methods have been used extensively in the neuroscience literature, they assume that there is unique information in each of a discrete set of frequency bands. The number of bands and frequency ranges used in these methods have been the subject of considerable controversy. Indeed Manning et al. (2009) have shown that broadband power often correlates more strongly with neuronal activity than does power at any narrow band. Also, nonparametric methods implicitly assume that the measured activity is observed independently during each observational epoch, and at each frequency, an assumption which is easily rejected in the data, which show strong temporal autocorrelation as well as correlations among frequency bands (von Stein and Sarnthein, 2000; Jensen and Colgin, 2007; Axmacher et al., 2010). Moreover, nonparametric methods are typically applied to EEG signals in a univariate fashion that neglects the spatial correlational structure. By simultaneously modeling the spatial and temporal structure in the data, parametric models confer greater statistical power so long as they are not poorly specified.
Parametric methods have been applied to various types of multivariate neural data including EEG (Hesse et al., 2003; Dhamala et al., 2008; Bastos et al., 2015), magnetoencephalography (MEG) (David et al., 2006), functional magnetic resonance imaging (FMRI) (Roebroeck et al., 2005; Goebel et al., 2003; David et al., 2008; Luo et al., 2013), and local field potentials (LFP) (Brovelli et al., 2004). These methods typically involve fitting vector autoregressive (VAR) models to multivariate neural data that are assumed to be stationary in a specific time interval of interest. The regression coefficient matrix derived from the VAR models can be used to study the flow of information between neuronal regions in the context of Granger causality (Gcausality). Neuroscientists have used Gaussian VAR models to study the effective connectivity (directed influence) between activated brain areas during cognitive and visuomotor tasks (Zhou et al., 2009; Deshpande et al., 2009; Graham et al., 2009; Roebroeck et al., 2005). Although VAR models and Gcausality methods have been argued to provide useful insights into the functional organization of the brain, their validity relies upon the assumptions of linearity and stationarity in mean and variance of the neural data. When one of these assumptions is violated, the conclusions drawn from a Gcausality analysis will be inconsistent and misleading (Seth et al., 2015). One of the most common violations by EEG signals is the assumption of variancestationarity (Wong et al., 2006). Therefore, in the present work, we adopt a stochastic volatility approach in which the nonstationary variance (also known as volatility) of the neural time series is assumed to follow a stochastic process. Such models have been extremely useful in the analyses of financial market data which, like neural data, exhibits high kurtosis (Heston, 1993; Bates, 1996; BarndorffNielsen and Shephard, 2002).
We propose a multivariate stochastic volatility (MSV) model with the aim of estimating the timevarying volatility of multivariate neural data and its spatial correlational structure. The MSV model assumes that the volatility series of iEEG signals follows a latentvariable vectorautoregressive process, and it allows for the lagged signals of different brain regions to influence each other by specifying a full persistent matrix (typically assumed to be diagonal) in the VAR process for volatility.
We employed a Bayesian approach to estimate the latent volatility series and the parameters of the MSV model using the forward filtering backward sampling and Metropolis Hastings algorithms. We validated the MSV model in a unique dataset comprising depthelectrode recordings from 96 neurosurgical patients. These patients volunteered to participate in a verbal recall memory task while they were undergoing clinical monitoring to localize the epileptogenic foci responsible for seizure onset. Our analyses focused on the subset of electrodes ($n=718$) implanted in medial temporal lobe (MTL) regions, including hippocampus, parahippocampal cortex, entorhinal cortex and perirhinal cortex. We chose to focus on these regions given their prominent role in the encoding and retrieval of episodic memories (Davachi et al., 2003; Kirwan and Stark, 2004; Kreiman et al., 2000; Sederberg et al., 2007).
We show that the MSV model, which allows for interactions between regions, provides a substantially superior fit to MTL recordings than univariate stochastic volatility (SV) models. The implied volatility in these models positively correlates with nonparametric estimates of spectral power, especially in the gamma frequency band.
We demonstrate the utility of our method for decoding cognitive states by using a penalized logistic regression classifier trained on the implied volatility data across MTL electrodes to predict which studied items will be subsequently recalled. We find that the MSVderived features outperform spectral features in decoding cognitive states, supporting the value of this modelbased timeseries analysis approach to the study of human cognition. Furthermore, using the MSV model to construct a directional MTL connectivity network, we find that significant bidirectional desynchronization between the perirhinal cortex and the hippocampus predicts successful memory encoding.
Multivariate Stochastic Volatility models for iEEG
Volatility of iEEG is stochastic
Previous studies have shown that variance of EEG recordings is timevarying (Wong et al., 2006; Galka et al., 2010). Klonowski et al. (2004) demonstrated a striking similarity between the timeseries of the Dow Jones index during economic recessions (big crashes) and the timeseries of iEEG during epileptic seizures. These timeseries typically possess 'big spikes' that are associated with abrupt changes in the variance of the measured signal. There are two main approaches for modeling the timevarying variance commonly known as volatility in the econometrics literature, both of which assume that volatility is a latent autoregressive process, that is the current volatility depends on its previous values. The first approach is the class of autoregressive conditional heteroscedastic (ARCH) models developed by Engle (1982) and the generalized autoregressive conditional heteroscedastic (GARCH) models extended by Bollerslev (1986). The ARCH/GARCH models assume that the current volatility is a deterministic function of the previous volatilities and the information up to the current time. The second approach is the class of stochastic volatility models proposed by Heston (1993), which assume that volatility is nondeterministic and follows a random process with a Gaussian error term. GARCHtype models have been popular in the empirical research community since the 1980’s due to their computational attractiveness. Inference of GARCH models is usually performed using maximum likelihood estimation methods (Engle and Bollerslev, 1986; Nelson, 1991). Until the mid 1990s, stochastic volatility models had not been widely used due to the difficulty in estimating the likelihood function, which involves integrating over a random latent process. With advances in computer technology, econometricians started to apply simulationbased techniques to estimate SV models (Kim et al., 1998; Chib et al., 2002; Jacquier et al., 1994). Despite their computational advantages, the GARCH model assumes a deterministic relationship between the current volatility and its previous information, making it slow to react to instantaneous changes in the system. In addition, Wang (2002) showed an asymptotic nonequivalance between the GARCH model and its diffusion limit if volatility is stochastic. Therefore, we employ a more general stochastic volatility approach to analyze iEEG signals.
To motivate the discussion of the MSV model, we first examine the distributional property of volatility. Previous econometrics literature has shown that volatilities of various financial timeseries can be well approximated by lognormal distributions such as the Standard and Poor 500 index (Cizeau et al., 1997), stock returns (Andersen, 2001a), and daily foreign exchange rates (Andersen et al., 2001b). The volatilities of these financial measures are found to be highly rightskewed and their logarithmic volatilities are approximately normal. Volatility of iEEG also exhibits this lognormality property. To demonstrate this property, we calculate and plot the density of the empirical variance of a sample iEEG timeseries using a rolling variance of window size 20 (Hull, 2009, chapter 17). Figure 1 demonstrates the distribution of the empirical volatility timeseries of the detrended (after removing autoregressive components) raw iEEG signals. The distribution is rightskewed and can be well approximated by a lognormal distribution. Due to this lognormality property, the SV approach typically models the logarithm of volatility instead of volatility itself. The next section describes the multivariate stochastic volatility model for iEEG data, which is a generalized version of the SV model that accounts for the interactions between recording brain locations.
The model
Stochastic volatility models belong to a wide class of nonlinear statespace models that have been extensively used in financial economics. There has been overwhelming evidence of nonstationarity in the variance of financial data (Black et al., 1972) and much effort has been made to model and understand the changes in volatility in order to forecast future returns, price derivatives, and study recessions, inflations and monetary policies (Engle and Patton, 2001; Cogley and Sargent, 2005; Blanchard and Simon, 2001). There is by now a large literature on stochastic volatility models and methods for estimating these models either by closedform solutions (Heston, 1993; Heston and Nandi, 2000) or by simulation (Harvey and Shephard, 1996; Kim et al., 1998; Omori et al., 2007). Under the stochastic volatility framework, the variance (or its monotonic transformation) of a time series is assumed to be a latent process, which is typically assumed to be autoregressive. Latent process models have also been widely applied to the neuroscience domain to model neuronal spiking activity using point processes (Macke et al., 2011; Smith and Brown, 2003; Eden et al., 2004), to study motor cortical activity using a statespace model (Wu et al., 2006), etc. These models have provided many insights into the mechanisms underlying cognitive processes. GARCHtype models have also been applied to EEG signals to study the transition into anesthesia in human, sleep stage transitions in sheep, and seizures in epileptic patients (Galka et al., 2010; Mohamadi et al., 2017). However, the applications of the GARCH models in the neuroscience literature remain on a smallscale, which focus on individual recording locations and neglect the connectivity among different regions of the brain. In this study, we provide a systematic way to study volatility of iEEG signals using a large iEEG dataset from the MTL region during a verbal memory task as a medium to illustrate how volatility and its connectivity network among MTL subregions can provide insights into the understanding of cognitive processes.
Following Harvey et al. (1994), we model the multivariate latent volatility process of the iEEG signals in the MTL region to follow a vector autoregressive model. The original model assumes that the coefficient matrix is diagonal, that is the past activity of one region does not have any influence on the others. In many financial applications, it is convenient to make this diagonality assumption to reduce the number of parameters of the MSV model, otherwise, a very large amount of data would be required to reliably estimate these parameters. We generalize the MSV model to allow for a full coefficient matrix in order to study the directional connections between different subregions in the MTL. This generalization is feasible due to hightemporalresolution neural time series collected using the intracranial electroencephalography. Let ${\mathbf{\text{y}}}_{t}=({y}_{1,t},\cdots ,{y}_{J,t})$ be a multivariate iEEG timeseries recordings at $J$ electrodes at time $t$. We model ${y}_{j,t}$, $1\le j\le J$, as follows:
and
where the error terms follow multivariate normal distributions: $\displaystyle {\mathit{\u03f5}}_{t}^{y}=({\u03f5}_{1,t}^{y},\cdots ,{\u03f5}_{J,t}^{y})\sim \mathcal{\mathcal{M}}\mathcal{\mathcal{V}}\mathcal{\mathcal{N}}(\mathbf{0},{I}_{J})$, ${\mathit{\u03f5}}_{t}^{x}=({\u03f5}_{1,t}^{x},\cdots ,{\u03f5}_{J,t}^{x})\sim \mathcal{\mathcal{M}}\mathcal{\mathcal{V}}\mathcal{\mathcal{N}}(\mathbf{0},\mathbf{\Sigma})$ denotes the identity matrix of rank $J$, and $\mathbf{\mathbf{\Sigma}}=diag({\sigma}_{1}^{2},\mathrm{\cdots},{\sigma}_{J}^{2})$ is assumed to be diagonal. That is, $\{{y}_{j,t}\}$ is a time series whose conditional logvariance (logvolatility), $\{{x}_{j,t}\}$, follows an AR(1) process that depends on its past value and the past values of other electrodes. The series {${y}_{1,t}\}{}_{t=1}{}^{T},\mathrm{\cdots},\{{y}_{J,t}\}{}_{t=1}{}^{T}$ are assumed to be conditionally independent given their logvolatility series ${\{{x}_{1,t}\}}_{t=1}^{T},\mathrm{\cdots},{\{{x}_{J,t}\}}_{t=1}^{T}$. The coefficient ${\beta}_{j,k}$ models how the past value of channel $k$ affects the current value of channel $j$ and ${\mu}_{k}$ is the unconditional average volatility at channel $k$. We can rewrite Equation 2 in a matrix form
where ${\text{\mathbf{x}}}_{t}=({x}_{1,t},\mathrm{\cdots},{x}_{J,t}),\mathit{\bm{\mu}}=({\mu}_{1},\mathrm{\cdots},{\mu}_{J})$, and $\mathit{\bm{\beta}}(j,k)={\beta}_{j,k}$. The vector error terms ${\mathit{\u03f5}}_{t}^{y}$ and ${\mathit{\u03f5}}_{t}^{x}$ are assumed to be independent. The parameters in the system above are assumed to be unknown and need to be estimated.
Following a Bayesian perspective, we assume that the parameters are not completely unknown, but they follow some prior distributions. Then, using the prior distributions and the information provided by the data, we can make inferences about the parameters from their posterior distributions.
Priors and estimation method
We specify prior distributions for the set of parameters $\mathit{\bm{\theta}}=(\mathit{\bm{\mu}},\mathit{\bm{\beta}},\mathbf{\mathbf{\Sigma}})$ of the MSV model. The mean vector $\mathit{\bm{\mu}}$ follows a flat multivariate normal distribution $\mathit{\bm{\mu}}\sim \mathcal{M}\mathcal{V}\mathcal{N}(0,1000{I}_{J})$. Each entry of the persistence matrix ${\beta}_{i,j}\in (1,1)$ is assumed to follow a beta distribution, $({\beta}_{i,j}+1)/2\sim Beta(20,1.5)$ (Kim et al., 1998). The beta prior distribution ensures that the entries of the persistent matrix are between −1 and 1, which guarantees the stationarity of the volatility process. For volatility of volatility, we utilize a flat gamma prior, ${\sigma}_{j}\sim \mathrm{\Gamma}(1/2,1/2\times 10)$ (Kastner and FrühwirthSchnatter, 2014), which is equivalent to $\pm \sqrt{{\sigma}_{j}^{2}}\sim N(0,10)$. We estimated the latent volatility processes and the parameters of the MSV model using a MetropoliswithinGibbs sampler (Kim et al., 1998; Omori et al., 2007; Kastner and FrühwirthSchnatter, 2014) (see Appendix 1 for derivation and 2 for discussion of parameter identification). Here, we choose hyperparameters which provide relatively flat prior distributions. In addition, due to a large amount of data (≈1 million time points per model for each subject) that was used to estimate the MSV model, the choice of the hyperparameters has little effect on the posterior estimates of the parameters in the MSV model.
Applications to verbalfree recall task
We analyzed the behavioral and electrophysiological data of 96 subjects implanted with MTL subdural and depth electrodes during a verbal free recall memory task (see Materials and methods section for details), a powerful paradigm for studying episodic memory. Subjects learned 25 lists of 12 unrelated words presented on a screen (encoding period) separated by 800–1200 ms interstimulus intervals (ISI), with each list followed by a short arithmetic distractor task to reduce recency effects (subjects more likely to recall words at the end of the list). During the retrieval period, subjects recalled as many words from the previously studied list as possible, in any order (Figure 2). In this paper, we focused our analyses on the MTL regions that have been implicated in episodic memory encoding (Squire and ZolaMorgan, 1991; Solomon et al., 2019; Long and Kahana, 2015). To assess a particular effect across subjects, we utilized the maximum a posteriori (MAP) estimate by taking the posterior mean of the variable of interest (whether it be the volatility time series ${\text{\mathbf{x}}}_{t}$ or the regression coefficient matrix $\mathit{\bm{\beta}}$) (Stephan et al., 2010).
Model comparison
To establish the validity of the MSV model, we compared its performance to that of univariate stochastic volatility models (equivalent to setting all the offdiagonal entries of the matrix $\mathit{\bm{\beta}}$ in Equation 2 to 0) in fitting iEEG data. We applied the MSV model to the multivariate neural data combined across encoding periods (regardless of whether the word items were later recalled) and SV models to datasets of individual electrodes with the assumption that the parameters of these models were shared across encoding events. We utilized the deviance information criterion (DIC) (Spiegelhalter et al., 2002; Gelman et al., 2014) considered to be a Bayesian analogue of the Akaike information criterion (AIC) to evaluate the performance of the models. The DIC consists of two components: the negative loglikelihood, $\overline{D}={\mathbb{E}}_{\mathit{\bm{\theta}},\text{\mathbf{x}}\mid \text{\mathbf{y}}}[2\mathrm{log}P(\text{\mathbf{y}}\mid \mathit{\bm{\theta}},\text{\mathbf{x}})]$, which measures the goodnessoffit of the model and the effective number of parameters, ${p}_{D}=\overline{D}D(\overline{\mathit{\bm{\theta}}},\overline{\text{\mathbf{x}}})={\mathbb{E}}_{\mathit{\bm{\theta}},\text{\mathbf{x}}\mid \text{\mathbf{y}}}[2\mathrm{log}P(\text{\mathbf{y}}\mid \mathit{\bm{\theta}},\text{\mathbf{x}})]+2\mathrm{log}P(\text{\mathbf{y}}\mid \overline{\mathit{\bm{\theta}}},\overline{\text{\mathbf{x}}})$, which measures the complexity of the model (see Appendix 5 for details on how to compute DIC for the MSV model). Where $\overline{\mathit{\bm{\theta}}}$ and $\overline{\text{\mathbf{x}}}$ denote the posterior means of the latent volatility series and the parameters of the MSV model. The DIC balances the tradeoff between model fit and model complexity. Models with smaller DICs are preferred. To account for the varying amount of data each subject had, we averaged the DIC by the number of events and electrodes. We found the MSV model to have a consistently lower DIC value than the SV model with a mean difference of 23 (±5.9 SEM). This indicates that the MSV model is approximately more than 150 times as probable as the SV models (Kass and Raftery, 1995), suggesting that the MSV model is a more appropriate model for iEEG data.
Relation to spectral power
We next analyzed the relation between volatility and spectral power (see Materials and methods) over a wide range of frequencies, from 3 to 180 Hz with 1 Hz steps). For each subject, we computed the correlation between volatility and spectral power for each encoding event and then averaged these correlations across all events. Since spectral powers of neighboring frequencies exhibit high correlations, we utilized a Gaussian regression model (Rasmussen, 2004) to estimate the correlation between volatility and spectral power as a function of frequency, which allows for a nonlinear relation. Figure 3 indicates that the correlation between volatility and spectral power is significantly positive across the spectrum and increasing in frequency. This illustrates the broadband nature of the volatility measure, but also suggests that volatility may more closely relate to previous neuroscientific findings observed for highfrequency as compared with lowfrequency activity. Having established that the MSV model outperforms the more traditional SV approach, and having shown that the implied volatility of the series reliably correlates with highfrequency neural activity, we next asked whether we can use the modelderived time series of volatility to predict subjects’ behavior in a memory task.
Classification of subsequent memory recall
Extensive previous work on the electrophysiological correlates of memory encoding has shown that spectral power, in both the lowfrequency (4–8 Hz) theta band and at frequencies about 40 Hz (socalled gamma activity), reliably predicts which studied words will be subsequently recalled or recognized (Sederberg et al., 2003). Here, we ask whether the implied volatility derived from the MSV model during word encoding can also reliably predict subsequent recall. To benchmark our MSV findings, we conducted parallel analyses of waveletderived spectral power at frequencies ranging between 3 and 180 Hz. To aggregate across MTL electrodes within each subject we applied an L2penalized logistic regression classifier using features extracted during the encoding period to predict subsequent memory performance (Ezzyat et al., 2017; Ezzyat et al., 2018). To estimate the generalization of the classifier, we utilized a nested crossvalidation procedure in which we trained the model on $N1$ sessions using the optimal penalty parameter selected via another inner crossvalidation procedure on the same training data (see Appendix 6 for details). We then tested the classifier on a holdout session collected at a different time. We computed the receiver operating characteristic (ROC) curve, relating true and false positives, as a function of the criterion used to assign regression output to response labels (see Appendix 6 for illustrations of ROC curves). We then use the AUC metric (area under the ROC curve) to characterize model performance. In order to perform a nested crossvalidation procedure, we focused on 42 subjects (out of 96 subjects with at least two sessions) with at least three sessions of recording. We find that MSVmodel implied volatility during item encoding reliably predicts subsequent recall, yielding an average AUC of 0.53 (95% CI, from 0.51 to 0.55). AUCs reliably exceeded chance levels in 72 percent of subjects (30 out of 42 subjects who contributed at least 3 sessions of data). Figure 4 compares these findings against results obtained using waveletderived power. Here we see that implied volatility does as well as, or better than, spectral measures at nearly all frequencies. In order to capture the correlation between spectral powers (thus their corresponding classifiers’ performances), we fit a Gaussian regression model to test the functional form of $\mathrm{\Delta}$AUC. We find that the $\mathrm{\Delta}$AUC function is significantly different from the 0 function (${\chi}_{11}^{2}=42$, P < ${10}^{5})$ (Benavoli and Mangili, 2015), which indicates that on average volatility performs significantly better than spectral power in predicting subsequent memory recall.
Directional connectivity analysis
Having established that volatility is predictive of subsequent memory recall, we now seek to identify directional connections between MTL subregions that are related to successful memory encoding. To investigate the intraMTL directional connectivity patterns that correlate with successful memory encoding, we utilize a subsequent memory effect (SME) paradigm in which we compare the MTL directional connectivity patterns (regression coefficient matrix $\mathit{\bm{\beta}}$) associated with recalled (R) word items to those associated with nonrecalled (NR) items. The SME paradigm has been widely used in the memory literature to study neural correlates (typically spectral power in a specific frequency band) that predict successful memory formation (Sederberg et al., 2003; Long et al., 2014; Burke et al., 2014). The intraMTL connectivity SME was constructed using the following procedure. First, we partitioned the word items into recalled and nonrecalled items offline. Using the MSV model, we constructed an intraMTL connectivity network for each memory outcome. We compared the distribution of the elements of these matrices across subjects. For the analysis, we considered four subregions of the MTL: hippocampus (Hipp), entorhinal cortex (EC), perirhinal cortex (PRC), and parahippocampal cortex (PHC). Each MTL subregion contains a different number of recording locations depending on the subject’s electrode coverage. We then computed the contrast between the two intraMTL networks corresponding to recalled and nonrecalled items for each ordered pair of subregions excluding the ones with fewer than 10 subjects contributing to the analysis. To compute the directional connectivity from region $I$ to region $J$, we took the average of the lagone 'influences’ that electrodes in region $I$ have on electrodes in region $J$, where $I$ denotes the number of electrodes in region $I$. We then computed the contrast between the two connectivity networks associated with recalled and nonrecalled items: ${\mathrm{\Delta}}_{I\to J}={C}_{I\to J}^{R}{C}_{I\to J}^{NR}$. Finally, we averaged the contrast for each ordered pair of MTL subregions across sessions within a subject. From now on, we refer to this contrast as simply connectivity network.
Figure 5 illustrates the intraMTL connectivity SME for the left and right hemispheres. Directed connections between the left hippocampus and the left PRC reliably decrease (falsediscoveryratecorrected) during successful memory encoding (${\mathrm{\Delta}}_{Hipp\to PRC}=0.04,{t}_{47}=3.49$, adj. P $<0.01$ and ${\mathrm{\Delta}}_{PRC\to Hipp}=0.06,{t}_{47}=2.66$, adj. P $<0.05)$. The difference between the directional connections between these two regions is not significant (${t}_{47}$=0.53, p=0.60). The decreases in the bidirectional connections within the left MTL are consistent with the findings in Solomon et al. (2017) which noted memoryrelated decreases in phase synchronization at high frequencies. We did not, however, find any other significant directional connections among the remaining regions (Figure 5, Appendix 7—tables 1 and 2).
Discussion
The ability to record electrophysiological signals from large numbers of brain recording sites has created a wealth of data on the neural basis of behavior and a pressing need for statistical methods suited to the properties of multivariate, neural, timeseries data. Because neural data strongly violate variancestationarity assumptions underlying standard approaches, such as Granger causality (Stokes and Purdon, 2017), researchers have generally eschewed these modelbased approaches and embraced nonparameter data analytic procedures. The multivariate stochastic volatility framework that we propose allows for nonstationary variance in the signals. This framework allows us to explicitly model the timevarying variance of neural signals. Similar stochastic volatility models have been used extensively in the financial economics literature to characterize a wide range of phenomena.
The MSV models proposed in this paper provide a new framework for studying multichannel neural data and relating them to cognition. Using a large MTL intracranical EEG dataset from 96 neurosurgical patients while performing a freerecall task, we found that volatility of iEEG timeseries is correlated with spectral power across the frequency spectrum and the correlation increases with frequency. To further test the ability of the MSV model to link iEEG recordings to behavior, we asked whether MTL volatility features during encoding can predict subsequent memory recall as well as spectral power features. Our findings indicate that volatility features significantly outperform the spectral features in decoding memory process in the human brain, suggesting that volatility can serve as a reliable measure for understanding cognitive processes.
A key strength of the MSV approach is its ability to identify directed interactions between brain regions without assuming variancestationarity of neural signals. We therefore used this approach to determine the directional connections between MTL subregions that correlate with successful memory encoding. Using the regression coefficient matrix of the multivariate volatility process, we found that periods of decreased connectivity in the volatility network among MTL subregions predicted successful learning. Specifically, we found that the hippocampus and the perirhinal cortex in the left hemisphere desynchronize (exerting less influence on one another) during successful learning, which is consistent with the latephase gamma decoupling noted in Fell et al. (2001). A more recent study by Solomon et al. (2019) also examined the association between intraMTL connectivity and successful memory formation using phasebased measures of connectivity. Solomon et al. (2019) noted intraMTL desynchronization at high frequencies during successful memory formation, aligning with the finding here that the volatility network tended to desynchronize. Furthermore, Solomon, et al. found broad increases in lowfrequency connectivity, which did not appear to be captured by our stochastic model. This suggests that volatility features reflect neural processes that are also captured by highfrequency phase information.
We further noted more statistically reliable changes in volatility networks in the left MTL compared to the right. This result is in line with a long history of neuroanatomical, electrophysiological, and imaging studies (e.g. Ojemann and Dodrill, 1985; Kelley et al., 1998) that found an association between verbal memory and the left MTL. It is possible that the verbal nature of our memory task specifically engaged processing in the left MTL, resulting in a lateralization of observed volatility phenomena.
Prior studies implicate the perirhinal cortex in judgement of familiarity and in recency discrimination system, while the hippocampus supports contextual binding (Eichenbaum et al., 2007; Diana et al., 2007; Hasselmo, 2005). These two systems play important roles in memory associative retrieval as suggested by animal studies (Brown and Aggleton, 2001), but it is still unclear how the hippocampus and PRC interact during memory processing. Fell et al. (2006) suggested that rhinalhippocampal coupling in the gamma range is associated with successful memory formation. Our results show no evidence for such a phenomenon, but rather agree with more recent studies demonstrating memoryrelated overall highfrequency desynchronization in the MTL (Solomon et al., 2019; Burke et al., 2015).
This paper presents the first major application of stochastic volatility models to neural timeseries data. The use of a multivariate modeling approach allows us to account for interactions between different subregions in the MTL and thus provides a better fit to the neural data than a univariate approach. Our MSV model fully captures how changes in neural data measured by volatility in one region influences changes in another region, providing insights into the complex dynamics of neural brain signals. We further demonstrated that volatility can be a promising biomarker due to its broadband nature by comparing its performance to one of spectral power in classifying subsequent memory. Finally, researchers can extend these models to broader classes of neural recordings, and exploit their statistical power to substantially increase our understanding of how behavior emerges from the complex interplay of neural activity across many brain regions.
Materials and methods
Participants
Ninety six patients with drugresistant epilepsy undergoing intracranial electroencephalographic monitoring were recruited in this study. Data were collected as part of a study of the effects of electrical stimulation on memoryrelated brain function at multiple medical centers. Surgery and iEEG monitoring were performed at the following centers: Thomas Jefferson University Hospital (Philadelphia, PA), Mayo Clinic (Rochester, MN), Hospital of the University of Pennsylvania (Philadelphia, PA), Emory University Hospital (Atlanta, GA), University of Texas Southwestern Medical Center (Dallas, TX), DartmouthHitchcock Medical Center (Lebanon, NH), Columbia University Medical Center (New York, NY) and the National Institutes of Health (Bethesda, MD). The research protocol was approved by the Institutional Review Board at each hospital and informed consent was obtained from each participant. Electrophysiological signals were collected from electrodes implanted subdurally on the cortical surface and within brain parenchyma. The neurosurgeons at each clinical site determined the placement of electrodes to best localize epileptogenic regions. Across the clinical sites, the following models of depth and grid electrodes (electrode diameter in parentheses) were used: PMT Depthalon (0.86 mm); Adtech Spencer RD (0.86 mm); Adtech Spencer SD (1.12 mm); Adtech BehnkeFried (1.28 mm); Adtech subdural and grids (2.3 mm). The dataset can be requested at http://memory.psych.upenn.edu/RAM_Public_Data.
Freerecall task
Request a detailed protocolEach subject participated in a delayed freerecall task in which they were instructed to study a list of words for later recall test. The task is comprised of three parts: encoding, delay, and retrieval. During encoding, the subjects were presented with a list of 12 words that were randomly selected from a pool of nouns (http://memory.psych.upenn.edu/WordPools). Each word presentation lasts for 1600 ms followed by a blank interstimulus interval (ISI) of 800 to 1200 ms. To mitigate the recency effect (recalling last items best) and the primacy effect (recalling first items better than the middle items), subjects were asked to perform a math distraction task immediately after the presentation of the last word. The math problems were of the form A+B+C = ?, where A,B,C were randomly selected digits. The delay math task lasted for 20 s, after which subjects were asked to recall as many words as possible from the recent list of words, in any order during the 30 s recall period. Subjects performed up to 25 lists per session of recording (300 words). Multiple sessions were recorded over the course of the patient’s hospital stay.
Electrophysiological recordings and data processing
Request a detailed protocoliEEG signals were recorded from subdural and depth electrodes at various sampling rates (500, 1000, or 1600 Hz) based on the the amplifier and the preference of the clinical team using one of the following EEG systems: DeltaMed XlTek (Natus), Grass Telefactor, and NihonKohden. We applied a 5 Hz bandstop fourth order Butterworth filter centered on 60 Hz to attenuate signal from electrical noise. We rereferenced the data using the common average of all electrodes in the MTL to eliminate potentially confounding largescale artifacts and noise. We used Morlet wavelet transform (wave number = 5) to compute power as a function of time for our iEEG signals. The frequencies were sample linearly from 3 to 180 Hz with 1 Hz increments. For each electrode and frequency, spectral power was logtransformed and then averaged over the encoding period. Within a session of recording, the spectral power was zscored using the distribution of power features across events.
To extract volatility feature, we applied the MSV model to the dataset constructed from the encoding events with an assumption that the parameters governing the dynamics of the volatility process does not change within a session of recording, that is the parameters of the MSV model are assumed to be shared across encoding events. Since we were only interested in the dynamics of the volatility time series of the brain signals, not the orignal time series themselves, we detrended the raw time series using vector autoregressive models of order $p$, where $p$ was selected based on the Akaike information criterion (AIC) to remove any autocorrelation in the raw signals and to make the time series more suited for an MSV application.
In the present manuscript, we used the common average reference (of MTL electrodes) to remove largescale noise from our MTL recordings. While the bipolar reference is frequently used for such analyses, due to its superior spatial selectivity, several factors limit its utility in this context. First, connectivity between a pair of adjacent bipolar electrodes is contaminated by signal common to their shared monopolar contact; as such, it is difficult to interpret connectivity between such pairs. In the setting of linear depth electrodes placed within the MTL, a substantial portion of the data between pairs of nearby MTL subregions would have to be excluded due to shared monopolar contacts. Second, bipolar rereferencing within the MTL produces the undesirable outcome that a bipolar midpoint 'virtual’ electrode could fall in a subregion where neither physical contact was placed, making observed connectivities difficult to interpret.
Anatomical localization
Request a detailed protocolThe MTL electrodes were anatomically localized using the following procedure. Hippocampal subfields and MTL cortices were automatically labeled in a preimplant 2 mm thick T2weighted MRI using the Automatic segmentation of hippocampal subfields (ASHS) multiatlas segmentation method (Yushkevich et al., 2015). A postimplant was coregistered with the MRI using Advanced Normalization Tools (Avants et al., 2008). MTL depth electrodes that were visible in the CT were then localized by a pair of neuroradiologists with expertise in MTL anatomy.
Statistical analyses
Request a detailed protocolTo assess an effect across subjects, we applied classical statistical tests on the maximum a posteriori (MAP) estimate of the parameter of interest . This approach has been used in many Bayesian applications to FMRI studies (Stephan et al., 2010) to test an effect across subjects. For analyses concerning frequencies, we applied Gaussian regression models (Rasmussen, 2004) to take the correlations among frequencies into account. We used the Matern (5/2) kernel function for all analyses that used Gaussian regression models, which enforces that the underlying functional form be at least twice differentiable. pvalues were FDRcorrected at $\alpha =0.05$ significance level when multiple tests were conducted.
Appendix 1
MCMC Algorithm for MSV Models
In this section, we derive an MCMC algorithm, which consists of MetropolisHastings steps within a Gibbs sampler, for estimating the latent volatility process and its parameters. As before, let $\mathbf{x}}_{t$ denote the latent multivariate volatility timeseries and $\theta =(\mu ,\beta ,\sum )$ the parameters of the MSV model. Following (Kim et al., 1998; Omori et al., 2007), we logtransform Equation 1
where $y}_{j,t}^{\ast$ is defined to be $\mathrm{l}\mathrm{o}\mathrm{g}({y}_{j,t}^{2}+c)$ with a fixed offset constant to $c={10}^{4}$ avoid values equal to 0. Equation 4 is linear but nonGaussian. To ameliorate the nonGaussianity problem, we approximate the logtransformed error term $\mathrm{log}\{({\u03f5}_{j,t}^{y}{)}^{2}\}\sim \mathrm{log}({\chi}_{1}^{2})$ by a mixture of 10 normal distributions as in Omori et al. (2007):
The values of ${p}_{k},{m}_{k}$ and ${v}_{k}$ are tabulated in Omori et al. (2007). As a result, we introduce a latent mixture component indicator variable, ${r}_{j,t}$, for channel $j$ at time $t$ such that $\mathrm{log}\{{({\u03f5}_{j,t}^{y})}^{2}\}\mid ({r}_{j,t}=k)\sim N({m}_{k},{v}_{k}^{2})$. The indicator variable is also estimated in the MCMC sampler. Given the mixture indicator ${\text{\mathbf{r}}}_{t}$ and the vector parameter $\mathit{\bm{\theta}}$, the latent volatility series ${\text{\mathbf{x}}}_{t}$ can be sampled using a forwardfiltering and backwardsampling (FFBS) procedure (West, 1996). The mixture indicator ${\text{\mathbf{r}}}_{t}$ can be sampled from a multinomial distribution
Finally, the vector parameter $\mathit{\bm{\theta}}$ can be sampled using an ancillaritysufficiency interweaving strategy (ASIS) (Yu and Meng, 2011; Kastner and FrühwirthSchnatter, 2014) which involves sampling the vector parameter $\mathit{\bm{\theta}}$ given the unstandardized volatility series ${x}_{j,t}$ via a MetropolisHasting step (noncentered step) and then sampling $\mathit{\bm{\theta}}$ again given the standardized volatility series ${\stackrel{~}{x}}_{j,t}={\displaystyle \frac{{x}_{j,t}{\mu}_{j}}{{\sigma}_{j}}}$ (centered step). Yu and Meng (2011) argued that by alternating between the noncentered and centered steps, we obtain a more efficient MCMC sampler that has a better mixing rate and converges faster. In addition, Kastner and FrühwirthSchnatter (2014) showed that the ASIS can accurately sample latent volatility timeseries that have low persistences, which is often the case for iEEG signals.
Appendix 2
Parameter Identification
Identification with Different Signaltonoise Ratios
To demonstrate that the Gibbs sampler can accurately estimate the latent volatility process and its associated parameters, we conducted a simulation study in which we generated $N=5$ time series of length $T=50,000$ according to Equations 1 and 2 with various signaltonoise ratios (SNR), which is controlled by varying the volatility of the volatility series, to mimic the typical length and the number of electrodes in our iEEG datasets. The SNR is calculated by taking the ratio of the average volatility of volatility across electrodes to the expected standard deviation of the noise term in Equation 4. We sampled 10,000 posterior draws and discarded the first 5000 draws as a burnin period to allow for convergence to the stationary distribution. Appendix 2—table 1 reports the posterior means of the parameters of the MSV model. Throughout the simulation, we use priors whose means are equal to the true values of the parameters. We observe that the Gibbs sampler can reliably estimate the parameters of the MSV model from datasets with various signaltonoise ratios. The identification of the parameters in the MSV model comes from the strength of our large iEEG dataset which typically has hundreds of thousands of data points per session, an amount of data that rarely exists in any financial application.
Appendix 3
How Much Data is Enough?
Here we provide an analysis of the amount of data needed to estimate the MSV model. The number of parameters of the MSV model is on the order of $O({J}^{2})$, where $J$ is the dimension of the multivariate time series. In our case, $J$ is just the number of recording locations. To justify this claim, we observe in Equation 2 that the variable reflecting the unconditional mean $\mathit{\bm{\mu}}$ has $J$ parameters, the persistent matrix $\mathit{\bm{\beta}}$ has ${J}^{2}$ parameters and the covariance matrix $\mathbf{\mathbf{\Sigma}}$ has $J$ parameters, totaling to $J+{J}^{2}+J={J}^{2}+2J=O({J}^{2})$ parameters. As the dimension of the timeseries increases, the amount of data required will need to increase quadratically for a good estimation of the model.
To assess the performance of various time lengths, we repeated the simulation study in Appendix .2 in which we simulated 30 datasets for each $T=k({J}^{2}+2J)$ for k ranging from 30 to 180 with an increment of 30. We used the logdeterminant error, which is defined to be
where ${\mathit{\bm{\beta}}}_{MSV}$ is the estimated connectivity matrix from the MSV model and ${\mathit{\bm{\beta}}}_{true}$ is the true simulated connectivity matrix. If ${\mathit{\bm{\beta}}}_{MSV}$ is close to ${\mathit{\bm{\beta}}}_{true}$, then the ${\mathit{\bm{\beta}}}_{MSV}^{1}\mathit{\bm{\beta}}$ is close to the identity matrix ${I}_{J}$, and therefore, the log determinant error will be very negative.
Appendix 3—figure 1 shows that as we increase the amount of data, the logdeterminant error becomes more negative. With $T=30({J}^{2}+2J)$, the average determinant error is approximately $8.5\times {10}^{9}$, which is very close . Therefore, we recommend the amount of data to be at least $30O({J}^{2})$ for the MSV model.
Appendix 3—figure 2 demonstrates that the MSV model can recover the latent volatility timeseries with $T=30({J}^{2}+2J)$, validating the algorithm that is used to estimate the paramters of the model.
Appendix 4
Model Fit Plots
We provide a visualization of the latent volatility series. Appendix 4—figure 1 illustrates the recordings from a hippocampal electrode during encoding of a list of 12 word items from a subject performing a verbal freerecall task at the University of Pennsylvania Hospital. The top panels show the detreneded iEEG series using AR(p) models. The bottom panels show the respective latent volatility series associated with the detrended signals. We observe that the latent volatility series capture the instantaneous variance of the original series.
Appendix 5
Model Selection Details
Deviance Information Criterion (DIC)
The deviance information criterion introduced in Spiegelhalter et al. (2002) is a Bayesian analogue of the Akaike information criterion (AIC), which is typically used to assess model performance. Given the observed data y and a parameter set $\mathit{\bm{\theta}}$ of a model, the AIC is defined to be
where $\widehat{\mathit{\bm{\theta}}}$ is the maximum likelihood estimate and $k$ is the number of parameters of the model. The AIC is just the negative loglikelihood (deviance) $D(\mathit{\bm{\theta}})=\mathrm{log}\{p(\text{\mathbf{y}}\mid \widehat{\mathit{\bm{\theta}}})$ penalized by the number of parameters used in the model. AIC has been shown to be asymptotically equivalent to leaveoneout crossvalidation (Stone, 1977). A small AIC indicates that the model fits the data well after accounting for the number of parameters. Under the MCMC framework, Spiegelhalter et al. (2002) proposed a Bayesian equivalence of Equation 7 such that the maximum likelihood estimate $\widehat{\mathit{\bm{\theta}}}$ is replaced by the posterior mean of $\stackrel{~}{\mathit{\bm{\theta}}}=E[\mathit{\bm{\theta}}\mid \text{\mathbf{y}}]$ and the number of parameters is replaced by the effective number of parameters ${p}_{D}={E}_{\mathit{\bm{\theta}}\mid \text{\mathbf{y}}}[2\mathrm{log}p(\text{\mathbf{y}}\mid \mathit{\bm{\theta}})]+2\mathrm{log}\{p(\text{\mathbf{y}}\mid \stackrel{~}{\mathit{\bm{\theta}}})\}$
The calculation of DIC can be readily done using the samples from the posterior distribution of $P(\mathit{\bm{\theta}}\mid \text{\mathbf{y}})$. Next we will derive the calculation of DIC for the MSV model.
DIC Calculation for MSV Model
The MCMC algorithm derived in Appendix 1 provides a way to sample from the posterior distribution of $\mathit{\bm{\theta}}=(\text{\mathbf{x}},\mathit{\bm{\mu}},\mathit{\bm{\beta}},\mathbf{\mathbf{\Sigma}})\mid \text{\mathbf{y}}$. We assume that there are $K$ (typically 1000) samples $\{{\mathit{\bm{\theta}}}^{(1)},\mathrm{\cdots}{\mathit{\bm{\theta}}}^{(K)}\}$ from the posterior distribution $P(\mathit{\bm{\theta}}\mid \text{\mathbf{y}})$. The posterior mean in Eqn. in eight can be estimated by $\stackrel{~}{\mathit{\bm{\theta}}}=(\stackrel{~}{\text{\mathbf{x}}},\stackrel{~}{\mathit{\bm{\mu}}},\stackrel{~}{\mathit{\bm{\beta}}},\stackrel{~}{\mathbf{\mathbf{\Sigma}}})={\displaystyle \frac{1}{K}}{\displaystyle \sum _{k=1}^{K}}{\mathit{\bm{\theta}}}^{(k)}$. The goodnessoffit term can be calculated as follows:
where $\varphi $ is the probability density function of the standard normal distribution. To calculate the effective number of parameters, we first need to compute the posterior mean of the loglikelihood function. This can be done using Monte Carlo integration
Appendix 6
Classification of Subsequent Memory Details
In this section, we explain the details of the penalized logistic regression classifier and the nested crossvalidation approach (Krstajic et al., 2014) that is used to evaluate the classifier’s predictive performance. We ask whether the implied volatility features across brain regions during the encoding period can reliably predict subsequent memory recall. Let ${\{{y}_{i}\}}_{i=1}^{N}$ be the binary outcomes (recalled versus nonrecalled) of $N$ wordtrials and ${\{{\text{\mathbf{X}}}_{i}\}}_{i=1}^{N}$ be their corresponding features of interest (volatility or spectral power). Here, we wish to distinguish ${y}_{i}$ and ${\text{\mathbf{X}}}_{i}$, which are the output and input variables of the classification model, from ${\{{\text{\mathbf{y}}}_{t}\}}_{t=1}^{T}$ and ${\{{\text{\mathbf{X}}}_{t}\}}_{t=1}^{T}$, which are iEEG and volatility timeseries respectively. The penalized L2 logistic classifier aims to find the optimal set of weights $\omega $ that minimizes the crossentropy between the outcome variable ${y}_{i}$ and the logistic transformation of ${\omega}^{T}{\text{\mathbf{X}}}_{i}$ while limiting the contributions of individual features by penalizing the L2norm of $\omega $
where ${\widehat{y}}_{i}={\displaystyle \frac{\mathrm{exp}({\omega}^{T}{\text{\mathbf{X}}}_{i})}{1+\mathrm{exp}({\omega}^{T}{\text{\mathbf{X}}}_{i})}}$ is the predicted probability of recall success, and $\lambda $ is the penalty parameter that controls the degree to which the classifier limits the contribution of individual features. A low value of $\lambda $ indicates that the classifier downgrades the contribution of the crossentropy in the loss function, and therefore penalizing the L2 norm of $\omega $ more severely and vice versa. In addition, minimizing the crossentropy loss between the empirical distribution (observed data) and distribution given by the model (predicted data) is equivalent to maximizing the loglikelihood function (Goodfellow et al., 2016). To test the ability of the classifier to generalize to new data, we employ a kfold nested crossvalidation approach in which we split the data naturally into different recording sessions (folds) illustrated in Appendix 6—figure 1. We hold out one session as the test data and train the classifier on the remaining session(s). The optimal hyperparameter $\lambda $ is optimized via another crossvalidation within the training set. The procedure is repeated until we obtain outofsample prediction for the entire dataset.
We assess the classifier’s performance using the area under the receiver operating curve (AUC), which relates the true positive rate and the false positive rate for various decision thresholds. An AUC of 1.0 implies that the classifier can perfectly discriminate between recalled and nonrecalled items and an AUC of 0.5 implies that the classifier is at chance. Appendix 6—figure 2 illustrates the classifier’s predictive performance for two subjects. The ROCs for these subjects sit above the 45degree line, indicating that the classifier can reliably discriminate subsequent memory recall above chance.
Appendix 7
Output Tables of Statistical Tests
This section reports the statistical tests for directional connection SME that contain at least 10 subjects. Appendix 7—tables 1 and 2 show the results of these tests for the left and right hemispheres.
References

1
The distribution of realized stock return volatilityJournal of Financial Economics 61:43–76.https://doi.org/10.1016/S0304405X(01)000551

2
The distribution of realized exchange rate volatilityJournal of the American Statistical Association 96:42–55.https://doi.org/10.1198/016214501750332965
 3
 4

5
Econometric analysis of realized volatility and its use in estimating stochastic volatility modelsJournal of the Royal Statistical Society: Series B 64:253–280.https://doi.org/10.1111/14679868.00336

6
EEG differences between eyesclosed and eyesopen resting conditionsClinical Neurophysiology 118:2765–2773.https://doi.org/10.1016/j.clinph.2007.07.028
 7

8
Jumps and stochastic volatility: exchange rate processes implicit in deutsche mark optionsReview of Financial Studies 9:69–107.https://doi.org/10.1093/rfs/9.1.69

9
Gaussian processes for bayesian hypothesis tests on regression functionsArtificial Intelligence and Statistics pp. 74–82.

10
The Capital Asset Pricing Model: Some Empirical TestsStanford University  Graduate School of Business.

11
The long and large decline in Us output volatilityBrookings Papers on Economic Activity 164:2001.https://doi.org/10.2139/ssrn.277356

12
Generalized autoregressive conditional heteroskedasticityJournal of Econometrics 31:307–327.https://doi.org/10.1016/03044076(86)900631
 13

14
Recognition memory: what are the roles of the perirhinal cortex and Hippocampus?Nature Reviews Neuroscience 2:51–61.https://doi.org/10.1038/35049064
 15

16
Human intracranial highfrequency activity during memory processing: neural oscillations or stochastic volatility?Current Opinion in Neurobiology 31:104–110.https://doi.org/10.1016/j.conb.2014.09.003
 17
 18

19
Human theta oscillations related to sensorimotor integration and spatial learningThe Journal of Neuroscience 23:4726–4736.https://doi.org/10.1523/JNEUROSCI.231104726.2003

20
Properties of slow oscillation during slowwave sleep and anesthesia in catsJournal of Neuroscience 31:14998–15008.https://doi.org/10.1523/JNEUROSCI.233911.2011

21
Markov chain monte carlo methods for stochastic volatility modelsJournal of Econometrics 108:281–316.https://doi.org/10.1016/S03044076(01)001373

22
Volatility distribution in the s&p500 stock indexPhysica A: Statistical Mechanics and Its Applications 4:441–445.https://doi.org/10.1016/S03784371(97)004172

23
Drifts and volatilities: monetary policies and outcomes in the post WWII USReview of Economic Dynamics 8:262–302.https://doi.org/10.1016/j.red.2004.10.009
 24
 25
 26

27
Multivariate granger causality analysis of fMRI dataHuman Brain Mapping 30:1361–1373.https://doi.org/10.1002/hbm.20606
 28

29
Imaging recollection and familiarity in the medial temporal lobe: a threecomponent modelTrends in Cognitive Sciences 11:379–386.https://doi.org/10.1016/j.tics.2007.08.001

30
Dynamic analysis of neural encoding by point process adaptive filteringNeural Computation 16:971–998.https://doi.org/10.1162/089976604773135069

31
The medial temporal lobe and recognition memoryAnnual Review of Neuroscience 30:123–152.https://doi.org/10.1146/annurev.neuro.30.051606.094328
 32

33
Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflationEconometrica: Journal of the Econometric Society, 1007, 10.2307/1912773.

34
Modelling the persistence of conditional variancesEconometric Reviews 5:1–50.https://doi.org/10.1080/07474938608800095

35
What good is a volatility model?Quantitative Finance 1:237–245.https://doi.org/10.1088/14697688/1/2/305
 36
 37

38
Human memory formation is accompanied by rhinalhippocampal coupling and decouplingNature Neuroscience 4:1259–1264.https://doi.org/10.1038/nn759
 39

40
Modeling Phase Transitions in the Brain27–52, Generalized statespace models for modeling nonstationary eeg timeseries, Modeling Phase Transitions in the Brain, Springer, 10.1007/9781441907967_2.

41
Understanding predictive information criteria for bayesian modelsStatistics and Computing 24:997–1016.https://doi.org/10.1007/s1122201394162
 42

43
Simultaneous EEG and fMRI of the alpha rhythmNeuroReport 13:2487–2492.https://doi.org/10.1097/0000175620021220000022
 44
 45

46
Multivariate stochastic variance modelsThe Review of Economic Studies 61:247–264.https://doi.org/10.2307/2297980

47
Estimation of an asymmetric stochastic volatility model for asset returnsJournal of Business & Economic Statistics 14:429–434.https://doi.org/10.1080/07350015.1996.10524672
 48

49
Cognitive functions of gammaband activity: memory match and utilizationTrends in Cognitive Sciences 8:347–355.https://doi.org/10.1016/j.tics.2004.06.006

50
The use of timevariant EEG Granger causality for inspecting directed interdependencies of neural assembliesJournal of Neuroscience Methods 124:27–44.https://doi.org/10.1016/S01650270(02)003667

51
A ClosedForm solution for options with stochastic volatility with applications to bond and currency optionsReview of Financial Studies 6:327–343.https://doi.org/10.1093/rfs/6.2.327

52
A ClosedForm GARCH option valuation modelReview of Financial Studies 13:585–625.https://doi.org/10.1093/rfs/13.3.585
 53

54
Direct brain recordings fuel advances in cognitive electrophysiologyTrends in Cognitive Sciences 14:162–171.https://doi.org/10.1016/j.tics.2010.01.005

55
Bayesian analysis of stochastic volatility modelsJournal of Business & Economic Statistics 12:371–389.https://doi.org/10.1080/07350015.1994.10524553

56
Crossfrequency coupling between neuronal oscillationsTrends in Cognitive Sciences 11:267–269.https://doi.org/10.1016/j.tics.2007.05.003

57
Theta returnsCurrent Opinion in Neurobiology 11:739–744.https://doi.org/10.1016/S09594388(01)002781

58
Bayes factorsJournal of the American Statistical Association 90:773–795.https://doi.org/10.1080/01621459.1995.10476572

59
Ancillaritysufficiency interweaving strategy (ASIS) for boosting MCMC estimation of stochastic volatility modelsComputational Statistics & Data Analysis 76:408–423.https://doi.org/10.1016/j.csda.2013.01.002
 60

61
Stochastic volatility: likelihood inference and comparison with ARCH modelsReview of Economic Studies 65:361–393.https://doi.org/10.1111/1467937X.00050
 62
 63

64
‘Epileptic seizures’ in economic organismPhysica A: Statistical Mechanics and Its Applications 342:701–707.https://doi.org/10.1016/j.physa.2004.05.045

65
Categoryspecific visual responses of single neurons in the human medial temporal lobeNature Neuroscience 3:946–953.https://doi.org/10.1038/78868
 66
 67

68
EEGcorrelated fMRI of human alpha activityNeuroImage 19:1463–1476.https://doi.org/10.1016/S10538119(03)002866
 69
 70
 71

72
Empirical models of spiking in neural populationsAdvances in Neural Information Processing Systems.
 73
 74

75
Arimagarch modeling for epileptic seizure prediction,IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE pp. 994–998.https://doi.org/10.1109/ICASSP.2017.7952305

76
Conditional heteroskedasticity in asset returns: a new approachEconometrica: Journal of the Econometric Society, 10.2307/2938260.
 77
 78

79
Stochastic volatility with leverage: fast and efficient likelihood inferenceJournal of Econometrics 140:425–449.https://doi.org/10.1016/j.jeconom.2006.07.008

80
Scikitlearn: machine learning in PythonJournal of Machine Learning Research 12:2825–2830.

81
Gating of human theta oscillations by a working memory taskThe Journal of Neuroscience 21:3175–3183.https://doi.org/10.1523/JNEUROSCI.210903175.2001

82
Advanced Lectures on Machine Learning63–71, Gaussian processes in machine learning, Advanced Lectures on Machine Learning, Springer, 10.1007/9783540286509_4.
 83

84
Theta and gamma oscillations during encoding predict subsequent recallThe Journal of Neuroscience 23:10809–10814.https://doi.org/10.1523/JNEUROSCI.233410809.2003
 85

86
Granger causality analysis in neuroscience and neuroimagingJournal of Neuroscience 35:3293–3297.https://doi.org/10.1523/JNEUROSCI.439914.2015

87
Estimating a StateSpace model from point process observationsNeural Computation 15:965–991.https://doi.org/10.1162/089976603765202622
 88
 89

90
Bayesian measures of model complexity and fitJournal of the Royal Statistical Society: Series B 64:583–639.https://doi.org/10.1111/14679868.00353

91
The medial temporal lobe memory systemScience 253:1380–1386.https://doi.org/10.1126/science.1896849

92
Ten simple rules for dynamic causal modelingNeuroImage 49:3099–3109.https://doi.org/10.1016/j.neuroimage.2009.11.015
 93

94
An asymptotic equivalence of choice of model by crossvalidation and akaike’s criterionJournal of the Royal Statistical Society. Series B 39:44–47.https://doi.org/10.1111/j.25176161.1977.tb01603.x

95
Different frequencies for different scales of cortical integration: from local gamma to long range alpha/theta synchronizationInternational Journal of Psychophysiology 38:301–313.https://doi.org/10.1016/S01678760(00)001720

96
Asymptotic nonequivalence of GARCH models and diffusionsThe Annals of Statistics 30:754–783.https://doi.org/10.1214/aos/1028674841
 97

98
Modelling nonstationary variance in EEG time series by state space GARCH modelComputers in Biology and Medicine 36:1327–1335.https://doi.org/10.1016/j.compbiomed.2005.10.001
 99

100
To center or not to center: that is not the Question—An Ancillarity–Sufficiency Interweaving Strategy (ASIS) for Boosting MCMC EfficiencyJournal of Computational and Graphical Statistics 20:531–570.https://doi.org/10.1198/jcgs.2011.203main
 101
 102
Decision letter

Michael J FrankSenior Editor; Brown University, United States

Michael BreakspearReviewing Editor; QIMR Berghofer Medical Research Institute, Australia

Pedro ValdésSosaReviewer
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Multivariate Stochastic Volatility Modeling of Neural Data" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Michael Frank as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Pedro ValdésSosa (Reviewer #3).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
The authors employ time varying analyses of the (multivariate) variance to a large data set of intracranial EEG. A hierarchical series of analyses is applied to capture an increasing complexity of distributed network activity. The authors establish that their approach can capture neuronal activity that predicts the behavioural performance of the experimental participants, an important proofofprinciple.
This is an important departure from classic (meanbased) autoregressive models of neurophysiological data. All three reviewers acknowledge the innovation of the approach and the quality of the presentation.
Essential revisions:
1) Better contextualization is requested by all three reviewers.
2) Reviewer 1 requests some further validation and I agree with this. You may choose to use eLife's figure supplements to do this.
3) Reviewer 1 also requests some elaboration of the theory and better presentation of the figures. These should be addressed.
4) I think the comparison with the power spectrum in the prediction section is informative in its current state (the authors already attended to aspects of this section in a presubmission process). Reviewer 1 has requested consideration of some alternative benchmarks here – these may be of interest but are not essential.
Reviewer #1:
The manuscript “Multivariate Stochastic Volatility Modeling of Neural Data” by Phan et al. addresses the problem of timevarying variance (volatility) in EEG data. It presents a multivariate model that describes the evolution of the variance over time. In this multivariate model, the variances of the channels evolve in a fashion described by an autoregressive model. The authors present a method to fit the model to data. They then (i) show that the multivariate model captures the data better than a univariate model, (ii) show that the volatility approach correlates with established methods based on the power spectrum, (iii) demonstrate that the volatility bears predictive power for the outcome of the experiment, and (iv) use the volatility coupling matrix to infer some changes in the effective connectivity between successfully memorized and not memorized items.
While their approach is intriguing in places, and the inference of multivariate couplings is indeed an open topic in the field, I have several considerations that limit my enthusiasm. I would request the authors to address these in detail before reconsidering the manuscript.
First, the authors introduce a novel method for estimating model parameters without showing many of its central properties. Most importantly, the datasensitivity of multivariate approaches is not explored here (although the authors do acknowledge it in the Introduction): how does the variability of their method depend on the amount of data? Appendix 1—table 1 does not address these problems. Furthermore, the three example studies are not a systematic proof of the consistency of the estimate. What happens, for example, if the largest eigenvalue of the coupling matrix is close to 1, i.e. close to a phase transition?
Second, in this light the authors should at least compare their method to established pairwise quantities (crosscorrelations, synchronization, information theoretic measures) and show that the observed decoupling is not equally visible in these wellestablished (and more dataefficient) methods.
Third, there is considerable literature on multivariate time series analysis, including volatility. For example, ARCH and GARCH models are rather frequently used. One reference that I am aware of (although this is just the tip of the iceberg) is Galka et al. (2010). The authors should relate their work to the established models, show how they differ and what benefit their novel approach yields.
Equation 1: I do not really understand why the model is constructed this way. I understand that the modeled variance enters exponentially, but do the authors have a justification for this approach? Second, the approach presented here only correlates the variance over time, while ϵjty is still drawn from a normal distribution. If I am not mistaken, y_{t} will hence look like Gaussian noise with changing amplitude. Does the data look like this? I would expect a more pronounced autocorrelation of the raw EEG signals. It might help if the authors showed (and clearly marked, see below about Figure 1) example data that they try to fit.
Figure 1: I find panel C inconsistent and incomprehensible in many respects. (i) Why are there 1 second intervals shown in the Countdown period? Why is this panel shown at all, given there are no data? (I guess this indicates that these data were not analyzed?) (ii) In the encoding period, why are there pauses between the words? These pauses are not mentioned in the description of the experiment. Are they included in the analysis or not? If not, why are the data still shown, given no data are shown for the Countdown period? Also, the label indicating the length of the intervals is rather hidden, given all other scales are at the top of the diagrams. Last, what does REM / Not REM mean? It can be guessed but is not explained anywhere. (iii) A distractor task following the encoding is mentioned in the text but not shown in the figure. (iv) The freerecall diagram is not explained at all. Was it analyzed (so far, the text says otherwise)? If not, why are data still shown, and what do the intervals indicate? Also there is no length scale for the white intervals. (v) The panel does not state what the plots show. The (very small, and unintuitively placed) label "Compute volatility" seems to indicate that it is the time series y_{t,i} of two channels, from which the volatility is then computed? (vi) I think readers could greatly benefit from an illustration of the workflow from y to theta and x, for example incorporating elements from Appendix 1—figure 1.
Subsection “Relation to Spectral Power”: This whole paragraph is imprecise in the details. (i) What is an "encoding event"? If it is each single word during any encoding period, this would imply a data length of 1.6 seconds each. How would this fit with the Introduction, which states that fitting a multivariate model is hard for limited data? (ii) The Materials and methods say that averaging was done per subject. So does Figure 2 show exemplary results for one subject? If so, this should be indicated, and the generalization to all subjects discussed.
Subsection “Classification of Subsequent Memory Recall”: I get that the authors want to show the predictive power of volatility. However, to make this point it would suffice to show that it is above chance level. I do not see the significance of comparing it to the predictive power of methods based on the power spectrum. Furthermore, why does the classifier only take into account each band of the power spectrum? Should a classification based on all bands simultaneously not be much more powerful?
Subsection “Directional Connectivity Analysis”: While I understand the approach by the authors, the results are displayed in a way that I find hard to comprehend. (i) The authors should indicate the inferred values for the coupling matrix and discuss its significance for the dynamics of the process (intrinsic timescales etc.). (ii) The authors first state that directed connections between left HC and PRC decrease, but then also say that the difference between directional connections between these areas is not significant. Do they refer to the difference between HC → PRC and PRC → HC? If yes, do the reported pvalues refer to the R or NR set of items?
The first three paragraphs of the Discussion are a mere introduction or repetition of results, and should hence be included in the Introduction or be omitted.
Discussion, fourth paragraph: Why is the decoupling only found in the left hemisphere? This striking difference is completely neglected by the authors. Although I am not too familiar with the cited study by Fell et al. (2006), it does not appear to confirm this difference between hemispheres. As this limitation to the left hemisphere will probably be very obvious to many readers, the authors should discuss it.
Discussion, last paragraph: The authors suggest that fellow researchers can use and extend their method. With this in mind, it would be very beneficial if the authors made their code publicly available.
Reviewer #2:
I believe that this paper succeeds in rolling out a new method for multivariate signal analysis for multichannel neural data. The paper presents an application of an approach to multiprocess (here multichannel) time series data that was first utilized in the analysis of financial time series. In this approach, the variance of the time series is not assumed to be constant across the periods of analysis. Methods such as Granger causality assume variance stationarity between brain regions in order to identify directed interactions between the regions. However, variance stationarity is often not observed in neural data. The method developed here (MSV: multivariate stochastic volatility) models multielectrode depth recordings from the medial temporal lobes (MTL) of 96 epilepsy patients participating in wordlist recall task. The key innovation is the use of an autoregressive framework to model the fluctuations in the variance of the time series across the encoding interval of the behavioral task and across the data channels. Through a fairly rigorous and detailed analysis of the variables of the statespace MSV model, including the regression coefficient matrix, the authors are able to infer how the signals in the multiple channels interact in time and across space. The authors are able to conclude that during successful encoding of words the connectivity between the hippocampus and the perirhinal cortex in the left hemisphere decreases from a resting condition where the former drives the latter more strongly than the latter drives the former. The authors are also able to predict which words from the lists will be successfully recalled by measuring the volatility in the neural data.
The methodology described here will find many applications in other studies using largescale electrophysiology and MEG. As a consequence, this paper will be of considerable interest to the large community of researchers investigating brain networks.
The authors use a limited number of equations to explain their methods of modeling and their analysis of the data and instead rely on the text do a lot of the expository work, which is fine given that text is clear and well written.
The authors do a good job of placing the MSV model in the context of nonparametric methods such as spectral analysis and provide a means to link the two together. Through this linkage the authors are able to show that the correlation between volatility and spectral power is nearly proportional to the analysis frequency in the spectrum. This suggests that neural activity in the gamma band is more volatile than at lower frequencies. This result is of interest because it suggests that the conflicting reports in the literature about the nature and function of gamma band activity in neural activity may in part be due to the application of nonparametric methods that require that the time series be stationary over the data periods of interest.
One improvement here would be if the authors also included the analysis of the neural activity during a time window other than the encoding phase. For example, if the countdown phase was also submitted to the same MSV modeling, the hippocampus and perirhinal cortex should show a different functional connectivity than what was reported for successful encoding of words. The authors should also comment on whether the MSV modeling can be applied to a network that includes 20 distinct sites, instead of just 4, as was done in the paper.
Reviewer #3:
This is excellent work and directs us away from looking at directional measure of influence based on pure autoregressive models (which focus on the conditional mean). Rather interest is focused on the conditional variance, which is shown in a clear way to have advantages over traditional Granger Causality measures. I am convinced by the evidence presented.
The paper would benefit by clarifying the subtypes of stochastic volatility models. There is some confusion in the field. The authors seem to be applying a type of GARCH model, sometimes proposed as distinct from SV and sometimes a subtype.
I do suggest that some prior work be cited properly. The Wong paper is cited in a way that suggests that it is relevant for emphasizing nonstationary when in reality is also a multivariate stochastic volatility model. The same group even attempted source localization in Galka, Yamashita and Ozaki (2004). This is all well covered in the book by Ozaki on time series modeling of neuroscience data.
Finally in the Discussion I would suggest a discussion with combined AR type models and those with stochastic volatility: see Mohamadi et al. (2017).
I also do not see if the code will be made publicly available.
https://doi.org/10.7554/eLife.42950.025Author response
Reviewer #1:
[…] While their approach is intriguing in places, and the inference of multivariate couplings is indeed an open topic in the field, I have several considerations that limit my enthusiasm. I would request the authors to address these in detail before reconsidering the manuscript.
First, the authors introduce a novel method for estimating model parameters without showing many of its central properties. Most importantly, the datasensitivity of multivariate approaches is not explored here (although the authors do acknowledge it in the Introduction): how does the variability of their method depend on the amount of data? Appendix 1—table 1 does not address these problems. Furthermore, the three example studies are not a systematic proof of the consistency of the estimate. What happens, for example, if the largest eigenvalue of the coupling matrix is close to 1, i.e. close to a phase transition?
Second, in this light the authors should at least compare their method to established pairwise quantities (crosscorrelations, synchronization, information theoretic measures) and show that the observed decoupling is not equally visible in these wellestablished (and more dataefficient) methods.
Third, there is considerable literature on multivariate time series analysis, including volatility. For example, ARCH and GARCH models are rather frequently used. One reference that I am aware of (although this is just the tip of the iceberg) is Galka et al. (2010). The authors should relate their work to the established models, show how they differ and what benefit their novel approach yields.
Equation 1: I do not really understand why the model is constructed this way. I understand that the modeled variance enters exponentially, but do the authors have a justification for this approach? Second, the approach presented here only correlates the variance over time, while ϵjty is still drawn from a normal distribution. If I am not mistaken, y_{t} will hence look like Gaussian noise with changing amplitude. Does the data look like this? I would expect a more pronounced autocorrelation of the raw EEG signals. It might help if the authors showed (and clearly marked, see below about Figure 1) example data that they try to fit.
We apologize for not providing enough details on the construction of the MSV model. The exponential term is actually the variance process, which is always positive. We model the logvariance using a vector autoregressive process since the variance is approximately lognormally distributed as demonstrated in the “Volatility of IEEG is Stochastic” section. In addition, we work with the detrended timeseries (after removing autoregressive components) instead of the original voltage series. The new Figure 1 demonstrates the process of going from the raw iEEG data to the volatility timeseries.
Figure 1: I find panel C inconsistent and incomprehensible in many respects. (i) Why are there 1 second intervals shown in the Countdown period? Why is this panel shown at all, given there are no data? (I guess this indicates that these data were not analyzed?). (ii) In the encoding period, why are there pauses between the words? These pauses are not mentioned in the description of the experiment. Are they included in the analysis or not? If not, why are the data still shown, given no data are shown for the Countdown period? Also, the label indicating the length of the intervals is rather hidden, given all other scales are at the top of the diagrams. Last, what does REM / Not REM mean? It can be guessed but is not explained anywhere. (iii) A distractor task following the encoding is mentioned in the text but not shown in the figure. (iv) The freerecall diagram is not explained at all. Was it analyzed (so far, the text says otherwise)? If not, why are data still shown, and what do the intervals indicate? Also there is no length scale for the white intervals. (v) The panel does not state what the plots show. The (very small, and unintuitively placed) label "Compute volatility" seems to indicate that it is the time series y_{t,i} of two channels, from which the volatility is then computed? (vi) I think readers could greatly benefit from an illustration of the workflow from y to theta and x, for example incorporating elements from Appendix 1—figure 1.
We replace Figure 1 in the previous manuscript with Figure 2 in the current manuscript. Figure 2 now includes the three relevant phases of the freerecall task, the recording locations, and the work flow of the MSV model.
Subsection “Relation to Spectral Power”: This whole paragraph is imprecise in the details. (i) What is an "encoding event"? If it is each single word during any encoding period, this would imply a data length of 1.6 seconds each. How would this fit with the Introduction, which states that fitting a multivariate model is hard for limited data? (ii) The Materials and methods say that averaging was done per subject. So does Figure 2 show exemplary results for one subject? If so, this should be indicated, and the generalization to all subjects discussed.
We apologize for not being clear. An encoding event is just a single word event. We apply the MSV model to the combined dataset of all encoding events with the assumption that these events share the same set of parameters which govern their volatility processes. This assumption is also stated in the Materials and methods section. Figure 2 in the previous manuscript now becomes Figure 3. In Figure 3, we clearly state that the confidence bands are constructed across subjects. Since each subject has multiple sessions of recordings, we run a separate MSV model for each session within a subject and the results are averaged across sessions.
Subsection “Classification of Subsequent Memory Recall”: I get that the authors want to show the predictive power of volatility. However, to make this point it would suffice to show that it is above chance level. I do not see the significance of comparing it to the predictive power of methods based on the power spectrum. Furthermore, why does the classifier only take into account each band of the power spectrum? Should a classification based on all bands simultaneously not be much more powerful?
We compare the performance of volatility to that of spectral power at a specific frequency in predicting subsequent memory recall to ensure that the two classifiers have the same number of features, thus making it a fair comparison. As you suggested, the classifier based on all bands yields a higher performance due to having a lot more information. What our manuscript suggests is that volatility contains more information than each of the individual frequencies.
Subsection “Directional Connectivity Analysis”: While I understand the approach by the authors, the results are displayed in a way that I find hard to comprehend. (i) The authors should indicate the inferred values for the coupling matrix and discuss its significance for the dynamics of the process (intrinsic timescales etc.). (ii) The authors first state that directed connections between left HC and PRC decrease, but then also say that the difference between directional connections between these areas is not significant. Do they refer to the difference between HC → PRC and PRC → HC? If yes, do the reported pvalues refer to the R or NR set of items?
(i) The inferred values reflect the lagone correlations among different locations in the brain, which are explained in the construction of the MSV model. (ii) Yes, the reviewer is correct. We compare the difference between HC→ PRCand PRC→ HC. However, each connection itself is a contrast between the recalleditem network and the nonrecalleditem network. We simply want to test to see whether the directional connections between these two regions are symmetric.
The first three paragraphs of the Discussion are a mere introduction or repetition of results, and should hence be included in the Introduction or be omitted.
We have fixed the first three paragraphs.
Discussion, fourth paragraph: Why is the decoupling only found in the left hemisphere? This striking difference is completely neglected by the authors. Although I am not too familiar with the cited study by Fell et al. (2006), it does not appear to confirm this difference between hemispheres. As this limitation to the left hemisphere will probably be very obvious to many readers, the authors should discuss it.
We include several citations in the Discussion section to explain this lateralization effect. In particular, the coupling only found in the left hemisphere is in line with neuroanatomical, electrophysiological, and imaging studies that found an association between verbal memory and the left MTL.
Discussion, last paragraph: The authors suggest that fellow researchers can use and extend their method. With this in mind, it would be very beneficial if the authors made their code publicly available.
The code has now been made available on GitHub.
Reviewer #2:
[…] The authors do a good job of placing the MSV model in the context of nonparametric methods such as spectral analysis and provide a means to link the two together. Through this linkage the authors are able to show that the correlation between volatility and spectral power is nearly proportional to the analysis frequency in the spectrum. This suggests that neural activity in the gamma band is more volatile than at lower frequencies. This result is of interest because it suggests that the conflicting reports in the literature about the nature and function of gamma band activity in neural activity may in part be due to the application of nonparametric methods that require that the time series be stationary over the data periods of interest.
One improvement here would be if the authors also included the analysis of the neural activity during a time window other than the encoding phase. For example, if the countdown phase was also submitted to the same MSV modeling, the hippocampus and perirhinal cortex should show a different functional connectivity than what was reported for successful encoding of words. The authors should also comment on whether the MSV modeling can be applied to a network that includes 20 distinct sites, instead of just 4, as was done in the paper.
As suggested by reviewer 2, we ran a connectivity analysis to compare the connections within the MTL region between the resting state period and the encoding period (including both recalled and nonrecalled items). However, we did not find any significant connections after correcting for multiple comparisons (see Author response image 1). This suggests that the contrast between recalled and nonrecalled items is perhaps stronger than the contrast between the resting state period and the encoding period.
Reviewer #3:
This is excellent work and directs us away from looking at directional measure of influence based on pure autoregressive models (which focus on the conditional mean). Rather interest is focused on the conditional variance, which is shown in a clear way to have advantages over traditional Granger Causality measures. I am convinced by the evidence presented.
The paper would benefit by clarifying the subtypes of stochastic volatility models. There is some confusion in the field. The authors seem to be applying a type of GARCH model, sometimes proposed as distinct from SV and sometimes a subtype.
I do suggest that some prior work be cited properly. The Wong paper is cited in a way that suggests that it is relevant for emphasizing nonstationary when in reality is also a multivariate stochastic volatility model. The same group even attempted source localization in Galka, Yamashita and Ozaki (2004). This is all well covered in the book by Ozaki on time series modeling of neuroscience data.
Finally in the Discussion I would suggest a discussion with combined AR type models and those with stochastic volatility: see Mohamadi et al. (2017).
I also do not see if the code will be made publicly available.
We add a literature review of volatility models in the “Volatility is Stochatic” section which includes a discussion of the GARCHtype models. In addition, we cite applications of GARCH models to EEG as suggested by reviewer 3. Furthermore, the MSV approach does incorporate ARtype model with stochastic volatility since we detrend the raw voltage timeseries first using an AR model and then apply the MSV model to the residual timeseries. The code of this study is now available on GitHub.
https://doi.org/10.7554/eLife.42950.026Article and author information
Author details
Funding
Defense Advanced Research Projects Agency (N660011424032)
 Michael J Kahana
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Blackrock Microsystems for providing neural recording and stimulation equipment. This work was supported by the DARPA Restoring Active Memory (RAM) program (Cooperative Agreement N660011424032). We owe a special thanks to the patients and their families for their participation and support of the study. The views, opinions, and/or findings contained in this material are those of the authors and should not be interpreted as representing the official views of the Department of Defense or the U.S. Government. MJK has started a company, Nia Therapeutics, LLC ('Nia'), intended to develop and commercialize brain stimulation therapies for memory restoration and has more than 5% equity interest in Nia. We thank Dr. James Kragel and Nicole Kratz for their thoughtful comments and inputs.
Ethics
Human subjects: Data were collected at the following centers: Thomas Jefferson University Hospital, Mayo Clinic, Hospital of the University of Pennsylvania, Emory University Hospital, University of Texas Southwestern Medical Center, DartmouthHitchcock Medical Center, Columbia University Medical Center, National Institutes of Health, and University of Washington Medical Center and collected and coordinated via the Data Coordinating Center (DCC) at the University of Pennsylvania. The research protocol for the Data Coordinating Center (DCC) was approved by the University of Pennsylvania IRB (protocol 820553) and informed consent was obtained from each participant.
Senior Editor
 Michael J Frank, Brown University, United States
Reviewing Editor
 Michael Breakspear, QIMR Berghofer Medical Research Institute, Australia
Reviewer
 Pedro ValdésSosa
Publication history
 Received: October 17, 2018
 Accepted: July 29, 2019
 Accepted Manuscript published: August 1, 2019 (version 1)
 Version of Record published: August 16, 2019 (version 2)
Copyright
© 2019, Phan et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 754
 Page views

 112
 Downloads

 1
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Download citations (links to download the citations from this article in formats compatible with various reference manager tools)
Open citations (links to open the citations from this article in various online reference manager services)
Further reading

 Computational and Systems Biology
 Physics of Living Systems

 Computational and Systems Biology
 Evolutionary Biology