Multivariate stochastic volatility modeling of neural data

Abstract
Introduction
Discussion
Materials and methods
Appendix 1
Appendix 2
Appendix 3
Appendix 4
Appendix 5
Appendix 6
Appendix 7
Data availability
References
Article and author information
Metrics

Abstract

Because multivariate autoregressive models have failed to adequately account for the complexity of neural signals, researchers have predominantly relied on non-parametric 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 latent-variable 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 high-frequency 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 perirhinal-hippocampal desynchronization in the MTL regions that is associated with successful memory encoding.

https://doi.org/10.7554/eLife.42950.001

Introduction

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 multi-electrode recordings, which yield a multivariate time series of voltage recorded at varying brain locations (denoted here as $𝐲_{t}$ ). Such signals may be measured non-invasively, 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 time-series 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 non-parametric spectral decomposition methods, as a means of uncovering features of neural activity that may correlate with behavior. A strength of these non-parametric methods is that they have enabled researchers to link fluctuations in iEEG signals to low-frequency neural oscillations observed during certain behavioral or cognitive states, such as slow-wave 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). High-frequency 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, non-parametric 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, non-parametric 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 (G-causality). 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 G-causality 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 G-causality analysis will be inconsistent and misleading (Seth et al., 2015). One of the most common violations by EEG signals is the assumption of variance-stationarity (Wong et al., 2006). Therefore, in the present work, we adopt a stochastic volatility approach in which the non-stationary 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; Barndorff-Nielsen and Shephard, 2002).

We propose a multivariate stochastic volatility (MSV) model with the aim of estimating the time-varying volatility of multivariate neural data and its spatial correlational structure. The MSV model assumes that the volatility series of iEEG signals follows a latent-variable vector-autoregressive 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 depth-electrode 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 non-parametric 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 MSV-derived features outperform spectral features in decoding cognitive states, supporting the value of this model-based time-series 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 time-varying (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 time-varying 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 non-deterministic and follows a random process with a Gaussian error term. GARCH-type 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 simulation-based 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 log-normal 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 right-skewed and their logarithmic volatilities are approximately normal. Volatility of iEEG also exhibits this log-normality 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 right-skewed and can be well approximated by a log-normal distribution. Due to this log-normality 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.

Figure 1

Download asset Open asset

Empirical characteristics of iEEG.

(A) Sample raw iEEG time series during a resting state (count down non-task) period for subject R1240T. (B) Detrended iEEG timeseries after removing autoregressive components. (C) Empirical variance timeseries calculated using a rolling-window of size 20. (D) Distribution of empirical variance with the blue curve showing the estimated empirical density using kernel density estimation and the red curve showing the best-fitting log-normal density to the data.

https://doi.org/10.7554/eLife.42950.002

The model

Stochastic volatility models belong to a wide class of non-linear state-space models that have been extensively used in financial economics. There has been overwhelming evidence of non-stationarity 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 closed-form 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 state-space model (Wu et al., 2006), etc. These models have provided many insights into the mechanisms underlying cognitive processes. GARCH-type 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 small-scale, 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 high-temporal-resolution neural time series collected using the intracranial electroencephalography. Let $y_{t} = (y_{1, t}, \dots, y_{J, t})$ be a multivariate iEEG time-series recordings at $J$ electrodes at time $t$ . We model $y_{j, t}$ , $1 \leq j \leq J$ , as follows:

y_{j, t} = \exp (\frac{x_{j, t}}{2}) ϵ_{j, t}^{y},

and

x_{j, t} - μ_{j} = \sum_{k = 1}^{J} β_{j, k} (x_{k, t - 1} - μ_{k}) + ϵ_{j, t}^{x},

where the error terms follow multivariate normal distributions: $ϵ_{t}^{y} = (ϵ_{1, t}^{y}, \dots, ϵ_{J, t}^{y}) \sim ℳ 𝒱 𝒩 (0, I_{J})$ , $ϵ_{t}^{x} = (ϵ_{1, t}^{x}, \dots, ϵ_{J, t}^{x}) \sim ℳ 𝒱 𝒩 (0, Σ)$ denotes the identity matrix of rank $J$ , and $𝚺 = d i a g (σ_{1}^{2}, \dots, σ_{J}^{2})$ is assumed to be diagonal. That is, ${y_{j, t}}$ is a time series whose conditional log-variance (log-volatility), ${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}, \dots, {y_{J, t}}_{t = 1}^{T}$ are assumed to be conditionally independent given their log-volatility series ${x_{1, t}}_{t = 1}^{T}, \dots, {x_{J, t}}_{t = 1}^{T}$ . The coefficient $β_{j, k}$ models how the past value of channel $k$ affects the current value of channel $j$ and $μ_{k}$ is the unconditional average volatility at channel $k$ . We can rewrite Equation 2 in a matrix form

𝐱_{t} - 𝝁 = 𝜷 (𝐱_{t - 1} - 𝝁) + ϵ_{t}^{x},

where $𝐱_{t} = (x_{1, t}, \dots, x_{J, t}), 𝝁 = (μ_{1}, \dots, μ_{J})$ , and $𝜷 (j, k) = β_{j, k}$ . The vector error terms $ϵ_{t}^{y}$ and $ϵ_{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 $𝜽 = (𝝁, 𝜷, 𝚺)$ of the MSV model. The mean vector $𝝁$ follows a flat multivariate normal distribution $𝝁 \sim ℳ 𝒱 𝒩 (0, 1000 I_{J})$ . Each entry of the persistence matrix $β_{i, j} \in (- 1, 1)$ is assumed to follow a beta distribution, $(β_{i, j} + 1) / 2 \sim B e t a (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, $σ_{j} \sim Γ (1 / 2, 1 / 2 \times 10)$ (Kastner and Frühwirth-Schnatter, 2014), which is equivalent to $\pm \sqrt{σ_{j}^{2}} \sim N (0, 10)$ . We estimated the latent volatility processes and the parameters of the MSV model using a Metropolis-within-Gibbs sampler (Kim et al., 1998; Omori et al., 2007; Kastner and Frühwirth-Schnatter, 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 verbal-free 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 Zola-Morgan, 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 $𝐱_{t}$ or the regression coefficient matrix $𝜷$ ) (Stephan et al., 2010).

Figure 2

Download asset Open asset

Task design and analysis.

(A) Subjects performed a verbal free-recall task which consists of three phases: (1) word encoding, (2) math distraction, and (3) retrieval. (B) 96 Participants were implanted with depth electrodes in the medial temporal lobe (MTL) with localized subregions: CA1, CA3, dentate gyrus (DG), subiculum (Sub), perirhinal cortex (PRC), entorhinal cortex (EC), or parahippocampal cortex (PHC). (C) To construct a directional connectivity network, we applied the MSV model to brain signals recorded from electrodes in the MTL during encoding. We analyzed the 1.6 s epochs during which words were presented on the screen. The network reflects directional lag-one correlations among the implied volatility timeseries recorded at various MTL subregions. (D) The upper row shows a sample of an individual patient’s raw voltage timeseries (blue) recorded from two electrodes during a word encoding period of 1.6 s, and the lower row shows their corresponding implied volatility timseries (gray) estimated using the MSV model.

https://doi.org/10.7554/eLife.42950.003

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 off-diagonal entries of the matrix $𝜷$ 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 log-likelihood, $\bar{D} = 𝔼_{𝜽, 𝐱 ∣ 𝐲} [- 2 \log P (𝐲 ∣ 𝜽, 𝐱)]$ , which measures the goodness-of-fit of the model and the effective number of parameters, $p_{D} = \bar{D} - D (\bar{𝜽}, \bar{𝐱}) = 𝔼_{𝜽, 𝐱 ∣ 𝐲} [- 2 \log P (𝐲 ∣ 𝜽, 𝐱)] + 2 \log P (𝐲 ∣ \bar{𝜽}, \bar{𝐱})$ , which measures the complexity of the model (see Appendix 5 for details on how to compute DIC for the MSV model). Where $\bar{𝜽}$ and $\bar{𝐱}$ denote the posterior means of the latent volatility series and the parameters of the MSV model. The DIC balances the trade-off 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 non-linear 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 high-frequency as compared with low-frequency 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 high-frequency neural activity, we next asked whether we can use the model-derived time series of volatility to predict subjects’ behavior in a memory task.

Figure 3

Download asset Open asset

Correlation between volatility and spectral power over a frequency range from 3 to 180 Hz.

We fit a Gaussian process model to estimate the functional form of the correlation function between volatility and spectral power (solid blue line). The 95% confidence bands were constructed from 96 subjects (dashed gray lines). The red line shows the null model. We observe a significantly positive correlation between volatility and spectral power, and the correlation increases with frequency.

https://doi.org/10.7554/eLife.42950.004

Classification of subsequent memory recall

Extensive previous work on the electrophysiological correlates of memory encoding has shown that spectral power, in both the low-frequency (4–8 Hz) theta band and at frequencies about 40 Hz (so-called 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 wavelet-derived spectral power at frequencies ranging between 3 and 180 Hz. To aggregate across MTL electrodes within each subject we applied an L2-penalized 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 cross-validation procedure in which we trained the model on $N - 1$ sessions using the optimal penalty parameter selected via another inner cross-validation procedure on the same training data (see Appendix 6 for details). We then tested the classifier on a hold-out 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 cross-validation 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 MSV-model 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 wavelet-derived 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 $Δ$ AUC. We find that the $Δ$ AUC function is significantly different from the 0 function ( $χ_{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.

Figure 4

Download asset Open asset

Classification of subsequent memory recall.

(A) Average AUC of the classifier trained on spectral power across 42 subjects with at least three sessions of recording (blue). The red line indicates the average AUC of the classifier trained on volatility. (B) $Δ$ AUC = AUC $_{v o l}$ $-$ AUC $_{p o w e r}$ as a function of frequency estimated by using a Gaussian regression model (dashed gray lines indicate 95% confidence bands). The red line shows the null model. We observe that the classifier trained on volatility performs at least as well as the one trained on spectral power across the frequency spectrum. We find that functional form of $Δ$ AUC is significantly different from the function ( $χ_{11}^{2} = 42$ , P < $10^{- 5}$ ) using a Gaussian process model, suggesting that the difference in performance between the volatility classifier and the spectral power classifier is significant.

https://doi.org/10.7554/eLife.42950.005

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 intra-MTL 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 $𝜷$ ) associated with recalled (R) word items to those associated with non-recalled (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 intra-MTL connectivity SME was constructed using the following procedure. First, we partitioned the word items into recalled and non-recalled items offline. Using the MSV model, we constructed an intra-MTL 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 intra-MTL networks corresponding to recalled and non-recalled 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 lag-one '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 non-recalled items: $Δ_{I \to J} = C_{I \to J}^{R} - C_{I \to J}^{N R}$ . 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 intra-MTL connectivity SME for the left and right hemispheres. Directed connections between the left hippocampus and the left PRC reliably decrease (false-discovery-rate-corrected) during successful memory encoding ( $Δ_{H i p p \to P R C} = - 0.04, t_{47} = - 3.49$ , adj. P $< 0.01$ and $Δ_{P R C \to H i p p} = - 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 bi-directional connections within the left MTL are consistent with the findings in Solomon et al. (2017) which noted memory-related 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).

Figure 5

Download asset Open asset

MTL directional connectivity network.

The MTL electrodes were divided into four subregions: hippocampus (Hipp.), parahippocampal cortex (PHC), entorhinal cortex (EC), and perirhinal cortex (PRC). The directional connectivity from region I to region J, $C_{I \to J} = \frac{1}{| I | | J |} \sum_{i \in I, j \in J} 𝜷_{i j}$ , was calculated by averaging the entries of the sub-matrix of the regression coefficient matrix $𝜷$ , whose rows and columns correspond to region I and J respectively. We computed the contrast between the directional connectivity of recalled and non-recalled events: $Δ_{I \to J} = C_{I \to J}^{R} - C_{I \to J}^{N R}$ for each subject. Solid lines show significant (FDR-corrected) connections between two regions and dashed lines show trending but insignificant connections. Red indicates positive changes and blue indicates negative changes. The directional connectivity from Hipp. to PRC is significant (adj. P < 0.01) and the reverse directional connectivity is also significant (adj. P < 0.05).

https://doi.org/10.7554/eLife.42950.006

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, time-series data. Because neural data strongly violate variance-stationarity assumptions underlying standard approaches, such as Granger causality (Stokes and Purdon, 2017), researchers have generally eschewed these model-based approaches and embraced non-parameter data analytic procedures. The multivariate stochastic volatility framework that we propose allows for non-stationary variance in the signals. This framework allows us to explicitly model the time-varying 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 multi-channel neural data and relating them to cognition. Using a large MTL intracranical EEG dataset from 96 neurosurgical patients while performing a free-recall 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 variance-stationarity 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 late-phase gamma decoupling noted in Fell et al. (2001). A more recent study by Solomon et al. (2019) also examined the association between intra-MTL connectivity and successful memory formation using phase-based measures of connectivity. Solomon et al. (2019) noted intra-MTL 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 low-frequency 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 high-frequency 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 rhinal-hippocampal 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 memory-related overall high-frequency 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 time-series 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

Key resources table

Reagent, type (species) or resource	Designation	Source or reference	Identifiers	Additional Information
Software and algorithm		Avants et al., 2008	http://picsl.upenn.edu/software/ants	advanced normalization tool
Software and algorithm		Yushkevich et al., 2015	https://www.nitrc.org/projects/ashs	ashs
Software and algorithm	sklearnc	Pedregosa et al., 2011	https://scikit-learn.org/stable/
Software and algorithm		This paper	http://memory.psych.upenn.edu/Electrophysiological_Data	custom processing scripts
Software and algorithm	PTSA	This paper	https://github.com/pennmem/ptsa_new	processing pipeline for reading in iEEG

Participants

Ninety six patients with drug-resistant 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 memory-related 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), Dartmouth-Hitchcock 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 Behnke-Fried (1.28 mm); Adtech subdural and grids (2.3 mm). The dataset can be requested at http://memory.psych.upenn.edu/RAM_Public_Data.

Free-recall task

Request a detailed protocol

Each subject participated in a delayed free-recall 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 inter-stimulus 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 protocol

iEEG 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 Nihon-Kohden. We applied a 5 Hz band-stop fourth order Butterworth filter centered on 60 Hz to attenuate signal from electrical noise. We re-referenced the data using the common average of all electrodes in the MTL to eliminate potentially confounding large-scale 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 log-transformed and then averaged over the encoding period. Within a session of recording, the spectral power was z-scored 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 large-scale 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 re-referencing 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 protocol

The MTL electrodes were anatomically localized using the following procedure. Hippocampal subfields and MTL cortices were automatically labeled in a pre-implant 2 mm thick T2-weighted MRI using the Automatic segmentation of hippocampal subfields (ASHS) multi-atlas segmentation method (Yushkevich et al., 2015). A post-implant was co-registered 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 protocol

To 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. p-values were FDR-corrected at $α = 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 Metropolis-Hastings steps within a Gibbs sampler, for estimating the latent volatility process and its parameters. As before, let $x_{t}$ denote the latent multivariate volatility time-series and $θ = (μ, β, \sum)$ the parameters of the MSV model. Following (Kim et al., 1998; Omori et al., 2007), we log-transform Equation 1

y_{j, t}^{*} = x_{j, t} + l o g {(ϵ_{j, t}^{y})^{2}},

where $y_{j, t}^{*}$ is defined to be $l o 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 non-Gaussian. To ameliorate the non-Gaussianity problem, we approximate the log-transformed error term $\log {(ϵ_{j, t}^{y})^{2}} \sim \log (χ_{1}^{2})$ by a mixture of 10 normal distributions as in Omori et al. (2007):

\log (χ_{1}^{2}) \sim \sum_{k = 1}^{10} P_{k} 𝒩 (m_{k}, v_{k}^{2}) .

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 $\log {{(ϵ_{j, t}^{y})}^{2}} ∣ (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 $𝐫_{t}$ and the vector parameter $𝜽$ , the latent volatility series $𝐱_{t}$ can be sampled using a forward-filtering and backward-sampling (FFBS) procedure (West, 1996). The mixture indicator $𝐫_{t}$ can be sampled from a multinomial distribution

P (r_{j, t} = k ∣ 𝐱_{t}, 𝜽) \propto P (r_{j, t} = k) \frac{1}{v_{k}} \exp {- \frac{{({\tilde{y}}_{j, t} - x_{j, t} - m_{k})}^{2}}{2 v_{k}^{2}}} .

Finally, the vector parameter $𝜽$ can be sampled using an ancillarity-sufficiency interweaving strategy (ASIS) (Yu and Meng, 2011; Kastner and Frühwirth-Schnatter, 2014) which involves sampling the vector parameter $𝜽$ given the unstandardized volatility series $x_{j, t}$ via a Metropolis-Hasting step (non-centered step) and then sampling $𝜽$ again given the standardized volatility series ${\tilde{x}}_{j, t} = \frac{x_{j, t} - μ_{j}}{σ_{j}}$ (centered step). Yu and Meng (2011) argued that by alternating between the non-centered and centered steps, we obtain a more efficient MCMC sampler that has a better mixing rate and converges faster. In addition, Kastner and Frühwirth-Schnatter (2014) showed that the ASIS can accurately sample latent volatility time-series that have low persistences, which is often the case for iEEG signals.

https://doi.org/10.7554/eLife.42950.008

Appendix 2

Parameter Identification

Identification with Different Signal-to-noise 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 signal-to-noise 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 burn-in 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 signal-to-noise 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 2—table 1

Parameter Recovery.

https://doi.org/10.7554/eLife.42950.010

		Truth							MSV
Dataset	SNR	Channel	$α$	$b_{1}$	$b_{2}$	$b_{3}$	$b_{4}$	$b_{5}$	$σ$	$α$	$b_{1}$	$b_{2}$	$b_{3}$	$b_{4}$	$b_{5}$	$σ$
1	0.16	1	3.39	0.85	0.20	0.00	0.00	0.00	0.17	3.39	0.86	0.19	0.00	−0.04	0.00	0.17
		2	3.60	0.00	0.88	−0.10	0.00	0.00	0.19	3.60	0.00	0.87	−0.10	−0.00	−0.01	0.19
		3	3.55	0.00	0.00	0.87	0.30	0.00	0.19	3.53	−0.00	−0.00	0.85	0.29	0.01	0.20
		4	3.51	0.00	0.00	0.00	0.71	0.00	0.12	3.51	0.01	−0.00	0.01	0.69	−0.03	0.14
		5	3.38	0.00	0.00	0.00	0.00	0.80	0.15	3.38	0.00	−0.01	−0.01	0.06	0.79	0.15
2	0.27	1	3.60	0.95	0.20	0.00	0.00	0.00	0.21	3.93	0.95	0.19	−0.00	0.00	0.00	0.24
		2	3.74	0.00	0.90	−0.10	0.00	0.00	0.25	3.82	−0.00	0.90	−0.10	0.00	−0.00	0.25
		3	3.89	0.00	0.00	0.93	0.30	0.00	0.29	3.81	−0.00	0.00	0.93	0.30	0.00	0.29
		4	3.05	0.00	0.00	0.00	0.91	0.00	0.28	3.04	−0.00	0.01	0.01	0.91	0.00	0.28
		5	3.96	0.00	0.00	0.00	0.00	0.95	0.22	3.97	−0.00	0.00	0.00	−0.00	0.94	0.22
3	0.42	1	3.25	0.52	0.20	0.00	0.00	0.00	0.42	3.25	0.51	0.20	0.01	−0.02	0.01	0.43
		2	3.64	0.00	0.57	−0.10	0.00	0.00	0.46	3.63	0.01	0.55	−0.14	0.02	0.02	0.46
		3	3.18	0.00	0.00	0.58	0.30	0.00	0.31	3.18	0.03	0.01	0.56	0.32	0.03	0.31
		4	3.44	0.00	0.00	0.00	0.65	0.00	0.41	3.45	0.00	−0.00	−0.00	0.64	0.00	0.41
		5	3.26	0.00	0.00	0.00	0.00	0.60	0.36	3.26	−0.02	−0.00	0.00	−0.01	0.60	0.36

We generated three datasets with different signal-to-noise ratios. The observed multivariate time-series $𝐲_{t}$ was simulated according to the data-generating process specified by the MSV model with pre-specified parameters (truth). We then applied the MSV model to the simulated series $𝐲_{t}$ to recover the parameters of the MSV model. In this simulation study, the non-zero off-diagonal entries of the matrix $𝜷$ were fixed across datasets. The diagonal elements of $𝜷$ were generated from a uniform distribution on $[0.7, 0.9]$ , $[0.9, 1.0]$ , and $[0.5, 0.7]$ respectively. The volatilities of volatility of the electrodes were generated from a uniform distribution on $[0.1, 0.2]$ , $[0.2, 0.3]$ , and $[0.3, 0.5]$ respectively.

https://doi.org/10.7554/eLife.42950.009

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 $𝝁$ has $J$ parameters, the persistent matrix $𝜷$ has $J^{2}$ parameters and the covariance matrix $𝚺$ has $J$ parameters, totaling to $J + J^{2} + J = J^{2} + 2 J = 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} + 2 J)$ for k ranging from 30 to 180 with an increment of 30. We used the log-determinant error, which is defined to be

log-determinant error = \log (| det (𝜷_{M S V}^{- 1} 𝜷_{t r u e} - I_{J}) |),

where $𝜷_{M S V}$ is the estimated connectivity matrix from the MSV model and $𝜷_{t r u e}$ is the true simulated connectivity matrix. If $𝜷_{M S V}$ is close to $𝜷_{t r u e}$ , then the $𝜷_{M S V}^{- 1} 𝜷$ 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 log-determinant error becomes more negative. With $T = 30 (J^{2} + 2 J)$ , 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 $30 O (J^{2})$ for the MSV model.

Appendix 3—figure 1

Download asset Open asset

Time length analysis.

For $k$ ranging from 30 to 180 with 30 increments and for each $T = k (J^{2} + 2 J)$ , we simulated 30 datasets according to the data-generating process specified by the MSV model. Then, we estimated the connectivity matrix $𝜷_{M S V}$ for each dataset. We assess the performance of the MSV model in estimating the true connectivity matrix $𝜷_{t r u e}$ using the log-determinant error metric, $\log (| det (𝜷_{M S V}^{- 1} 𝜷_{t r u e} - I_{J}) |)$ . The figure shows the average performance at each time length with the standard error bars.

https://doi.org/10.7554/eLife.42950.012

Appendix 3—figure 2 demonstrates that the MSV model can recover the latent volatility timeseries with $T = 30 (J^{2} + 2 J)$ , validating the algorithm that is used to estimate the paramters of the model.

Appendix 3—figure 2

Download asset Open asset

Volatility Timeseries Recovery.

The red lines show the true simulated log-volatility timeseries for the first 1000 time points for two channels . The blue timeseries show the estimated log-volatility time series using the MSV model and the gray timeseries are the 95% posterior confidence intervals. The figure demonstrates that the MSV model can estimate the latent log-volatility timeseries well.

https://doi.org/10.7554/eLife.42950.013

https://doi.org/10.7554/eLife.42950.011

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 free-recall 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 4—figure 1

Download asset Open asset

Model fit plots for a hippocampal electrode.

The upper panels show the detrended iEEG signals using an AR(p) model for encoding periods of a list of words. The lower panels show the associated estimated latent volatility processes. The red lines indicate the average volatility during the encoding period. Red words are later recalled and blue words are not recalled.

https://doi.org/10.7554/eLife.42950.015

https://doi.org/10.7554/eLife.42950.014

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 $𝜽$ of a model, the AIC is defined to be

AIC = \underset{goodness of fit}{\underset{⏟}{- 2 \log {p (𝐲 ∣ \hat{𝜽})}}} + \underset{number of parameters}{\underset{⏟}{2 k}},

where $\hat{𝜽}$ is the maximum likelihood estimate and $k$ is the number of parameters of the model. The AIC is just the negative log-likelihood (deviance) $D (𝜽) = \log {p (𝐲 ∣ \hat{𝜽})$ penalized by the number of parameters used in the model. AIC has been shown to be asymptotically equivalent to leave-one-out cross-validation (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 $\hat{𝜽}$ is replaced by the posterior mean of $\tilde{𝜽} = E [𝜽 ∣ 𝐲]$ and the number of parameters is replaced by the effective number of parameters $p_{D} = E_{𝜽 ∣ 𝐲} [- 2 \log p (𝐲 ∣ 𝜽)] + 2 \log {p (𝐲 ∣ \tilde{𝜽})}$

DIC = \underset{goodness of fit}{\underset{⏟}{- 2 \log {p (𝐲 ∣ \tilde{𝜽})}}} + \underset{effective number of parameters}{\underset{⏟}{2 p_{D}}} .

The calculation of DIC can be readily done using the samples from the posterior distribution of $P (𝜽 ∣ 𝐲)$ . 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 $𝜽 = (𝐱, 𝝁, 𝜷, 𝚺) ∣ 𝐲$ . We assume that there are $K$ (typically 1000) samples ${𝜽^{(1)}, \dots 𝜽^{(K)}}$ from the posterior distribution $P (𝜽 ∣ 𝐲)$ . The posterior mean in Eqn. in eight can be estimated by $\tilde{𝜽} = (\tilde{𝐱}, \tilde{𝝁}, \tilde{𝜷}, \tilde{𝚺}) = \frac{1}{K} \sum_{k = 1}^{K} 𝜽^{(k)}$ . The goodness-of-fit term can be calculated as follows:

- 2 \log p (y ∣ \tilde{θ}) = - 2 \log p (y ∣ \tilde{x}, \tilde{μ}, \tilde{β}, \tilde{Σ})

= - 2 \sum_{t = 1}^{T} \log p (y_{t} ∣ {\tilde{x}}_{t})

= - 2 \sum_{t = 1}^{T} \log ϕ (y_{t} \exp (- \tilde{x_{t}} / 2)),

where $ϕ$ 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 log-likelihood function. This can be done using Monte Carlo integration

E_{θ ∣ y} [- 2 \log p (y ∣ θ)] = \frac{- 2}{K} \sum_{k = 1}^{K} \log p (y ∣ θ^{(k})

= \frac{- 2}{K} \sum_{k = 1}^{K} \sum_{t = 1}^{T} \log ϕ (y_{t} \exp (- x_{t}^{(k)} / 2)) .

https://doi.org/10.7554/eLife.42950.016

Appendix 6

Classification of Subsequent Memory Details

In this section, we explain the details of the penalized logistic regression classifier and the nested cross-validation 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 non-recalled) of $N$ word-trials and ${𝐗_{i}}_{i = 1}^{N}$ be their corresponding features of interest (volatility or spectral power). Here, we wish to distinguish $y_{i}$ and $𝐗_{i}$ , which are the output and input variables of the classification model, from ${𝐲_{t}}_{t = 1}^{T}$ and ${𝐗_{t}}_{t = 1}^{T}$ , which are iEEG and volatility timeseries respectively. The penalized L-2 logistic classifier aims to find the optimal set of weights $ω$ that minimizes the cross-entropy between the outcome variable $y_{i}$ and the logistic transformation of $ω^{T} 𝐗_{i}$ while limiting the contributions of individual features by penalizing the L2-norm of $ω$

l (ω) = - λ \underset{cross-entropy}{\underset{⏟}{\sum_{i = 1}^{N} y_{i} \log \hat{y_{i}} + (1 - y_{i}) \log (1 - {\hat{y}}_{i})}} + \underset{penalty}{\underset{⏟}{0.5 {|| ω ||}_{2}^{2}}},

where ${\hat{y}}_{i} = \frac{\exp (ω^{T} 𝐗_{i})}{1 + \exp (ω^{T} 𝐗_{i})}$ is the predicted probability of recall success, and $λ$ is the penalty parameter that controls the degree to which the classifier limits the contribution of individual features. A low value of $λ$ indicates that the classifier downgrades the contribution of the cross-entropy in the loss function, and therefore penalizing the L-2 norm of $ω$ more severely and vice versa. In addition, minimizing the cross-entropy loss between the empirical distribution (observed data) and distribution given by the model (predicted data) is equivalent to maximizing the log-likelihood function (Goodfellow et al., 2016). To test the ability of the classifier to generalize to new data, we employ a k-fold nested cross-validation 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 $λ$ is optimized via another cross-validation within the training set. The procedure is repeated until we obtain out-of-sample prediction for the entire dataset.

Appendix 6—figure 1

Download asset Open asset

Cross-validation scheme.

For each subject, we train our classifier on all (yellow) but one session and test the performance on the hold-out session (green) and repeat the procedure for each session.

https://doi.org/10.7554/eLife.42950.018

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 non-recalled 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 45-degree line, indicating that the classifier can reliably discriminate subsequent memory recall above chance.

Appendix 6—figure 2

Download asset Open asset

Sample ROC curves (gray).

The black dashed lines indicate ROC curves of at-chance classifiers. The area under the ROC curve (AUC) is a measure of the performance of the classifier across a spectrum of decision thresholds. An AUC of 1 indicates that the classifier can classify recalled items perfectly. An AUC of 0.5 indicates that the classifier is at chance, that is as good as a random coin flip.

https://doi.org/10.7554/eLife.42950.019

https://doi.org/10.7554/eLife.42950.017

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.

Appendix 7—table 1

intra-MTL directional connectivity in the left hemisphere.

https://doi.org/10.7554/eLife.42950.021

Region I $\to$ region J	Mean $Δ_{I \to J}$	Se	T	N	P	adj. P
Hipp $\to$ PRC	−0.044	0.013	−3.494	48	0.001	0.009**
PRC $\to$ Hipp	−0.060	0.022	−2.667	48	0.011	0.045*
Hipp $\to$ EC	0.010	0.032	0.312	14	0.768	0.953
EC $\to$ Hipp	−0.077	0.067	−1.158	14	0.285	0.569
PHC $\to$ PRC	−0.006	0.040	−0.146	16	0.889	0.953
PRC $\to$ PHC	−0.037	0.028	−1.348	16	0.212	0.564
PRC $\to$ EC	0.001	0.022	0.061	21	0.953	0.953
EC $\to$ PRC	0.005	0.030	0.168	21	0.872	0.953

Appendix 7—table 2

intra-MTL directional connectivity in the right hemisphere.

https://doi.org/10.7554/eLife.42950.022

Region I $\to$ region J	Mean $Δ_{I \to J}$	Se	T	N	P	adj. P
Hipp $\to$ PRC	−0.010	0.020	−0.471	40	0.645	0.838
PRC $\to$ Hipp	−0.016	0.029	−0.575	40	0.574	0.838
Hipp $\to$ EC	0.011	0.030	0.361	14	0.733	0.838
EC $\to$ Hipp	0.054	0.047	1.144	14	0.290	0.838
Hipp $\to$ PHC	−0.027	0.032	−0.837	15	0.432	0.838
PHC $\to$ Hipp	0.044	0.035	1.232	15	0.254	0.838
PRC $\to$ EC	0.020	0.052	0.378	14	0.722	0.838
EC $\to$ PRC	−0.002	0.081	−0.020	14	0.985	0.985

$*, :$ P $< 0.05$ , $*, *, :$ P $< 0.01$ , adjusted $p$ -values were calculated using the Benjamini-Hochberg procedure.

https://doi.org/10.7554/eLife.42950.020

Data availability

The iEEG dataset collected from epileptic patients in this paper is available and, to protect patients' confidentiality, can be requested at http://memory.psych.upenn.edu/RAM_Public_Data. The cmlreaders repository for reading in the data is at https://github.com/pennmem/. The main script for the paper is available at https://github.com/tungphan87/MSV_EEG (copy archived at https://github.com/elifesciences-publications/MSV_EEG).

References

1. Andersen T
(2001a) The distribution of realized stock return volatility
Journal of Financial Economics 61:43–76.

https://doi.org/10.1016/S0304-405X(01)00055-1
- Google Scholar
(2001b) The distribution of realized exchange rate volatility
Journal of the American Statistical Association 96:42–55.

https://doi.org/10.1198/016214501750332965
- Google Scholar
(2008) Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain
Medical Image Analysis 12:26–41.

https://doi.org/10.1016/j.media.2007.06.004
- PubMed
- Google Scholar
1. Axmacher N
2. Henseler MM
3. Jensen O
4. Weinreich I
5. Elger CE
6. Fell J
(2010) Cross-frequency coupling supports multi-item working memory in the human Hippocampus
PNAS 107:3228–3233.

https://doi.org/10.1073/pnas.0911531107
- PubMed
- Google Scholar
1. Barndorff-Nielsen OE
2. Shephard N
(2002) Econometric analysis of realized volatility and its use in estimating stochastic volatility models
Journal of the Royal Statistical Society: Series B 64:253–280.

https://doi.org/10.1111/1467-9868.00336
- Google Scholar
(2007) EEG differences between eyes-closed and eyes-open resting conditions
Clinical Neurophysiology 118:2765–2773.

https://doi.org/10.1016/j.clinph.2007.07.028
- PubMed
- Google Scholar
1. Bastos AM
2. Vezoli J
3. Bosman CA
4. Schoffelen JM
5. Oostenveld R
6. Dowdall JR
7. De Weerd P
8. Kennedy H
9. Fries P
(2015) Visual Areas exert feedforward and feedback influences through distinct frequency channels
Neuron 85:390–401.

https://doi.org/10.1016/j.neuron.2014.12.018
- PubMed
- Google Scholar
1. Bates DS
(1996) Jumps and stochastic volatility: exchange rate processes implicit in deutsche mark options
Review of Financial Studies 9:69–107.

https://doi.org/10.1093/rfs/9.1.69
- Google Scholar
1. Benavoli A
2. Mangili F
(2015)
Gaussian processes for bayesian hypothesis tests on regression functions

Artificial Intelligence and Statistics pp. 74–82.
- Google Scholar
Book
(1972)
The Capital Asset Pricing Model: Some Empirical Tests

Stanford University - Graduate School of Business.
- Google Scholar
1. Blanchard O
2. Simon J
(2001) The long and large decline in Us output volatility
Brookings Papers on Economic Activity 164:2001.

https://doi.org/10.2139/ssrn.277356
- Google Scholar
1. Bollerslev T
(1986) Generalized autoregressive conditional heteroskedasticity
Journal of Econometrics 31:307–327.

https://doi.org/10.1016/0304-4076(86)90063-1
- Google Scholar
1. Brovelli A
2. Ding M
3. Ledberg A
4. Chen Y
5. Nakamura R
6. Bressler SL
(2004) Beta oscillations in a large-scale sensorimotor cortical network: directional influences revealed by granger causality
PNAS 101:9849–9854.

https://doi.org/10.1073/pnas.0308538101
- PubMed
- Google Scholar
1. Brown MW
2. Aggleton JP
(2001) Recognition memory: what are the roles of the perirhinal cortex and Hippocampus?
Nature Reviews Neuroscience 2:51–61.

https://doi.org/10.1038/35049064
- Google Scholar
1. Burke JF
2. Long NM
3. Zaghloul KA
4. Sharan AD
5. Sperling MR
6. Kahana MJ
(2014) Human intracranial high-frequency activity maps episodic memory formation in space and time
NeuroImage 85 Pt 2:834–843.

https://doi.org/10.1016/j.neuroimage.2013.06.067
- PubMed
- Google Scholar
(2015) Human intracranial high-frequency 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
- PubMed
- Google Scholar
(2007) Remembering the past and imagining the future: a neural model of spatial memory and imagery
Psychological Review 114:340–375.

https://doi.org/10.1037/0033-295X.114.2.340
- Google Scholar
1. Canolty RT
2. Edwards E
3. Dalal SS
4. Soltani M
5. Nagarajan SS
6. Kirsch HE
7. Berger MS
8. Barbaro NM
9. Knight RT
(2006) High gamma power is phase-locked to theta oscillations in human neocortex
Science 313:1626–1628.

https://doi.org/10.1126/science.1128115
- PubMed
- Google Scholar
(2003) Human theta oscillations related to sensorimotor integration and spatial learning
The Journal of Neuroscience 23:4726–4736.

https://doi.org/10.1523/JNEUROSCI.23-11-04726.2003
- PubMed
- Google Scholar
(2011) Properties of slow oscillation during slow-wave sleep and anesthesia in cats
Journal of Neuroscience 31:14998–15008.

https://doi.org/10.1523/JNEUROSCI.2339-11.2011
- PubMed
- Google Scholar
(2002) Markov chain monte carlo methods for stochastic volatility models
Journal of Econometrics 108:281–316.

https://doi.org/10.1016/S0304-4076(01)00137-3
- Google Scholar
1. Cizeau P
2. Liu Y
3. Meyer M
4. Peng C-K
5. Stanley HE
(1997) Volatility distribution in the s&p500 stock index
Physica A: Statistical Mechanics and Its Applications 4:441–445.

https://doi.org/10.1016/S0378-4371(97)00417-2
- Google Scholar
1. Cogley T
2. Sargent TJ
(2005) Drifts and volatilities: monetary policies and outcomes in the post WWII US
Review of Economic Dynamics 8:262–302.

https://doi.org/10.1016/j.red.2004.10.009
- Google Scholar
(2003) Multiple routes to memory: distinct medial temporal lobe processes build item and source memories
PNAS 100:2157–2162.

https://doi.org/10.1073/pnas.0337195100
- PubMed
- Google Scholar
1. David O
2. Kiebel SJ
3. Harrison LM
4. Mattout J
5. Kilner JM
6. Friston KJ
(2006) Dynamic causal modeling of evoked responses in eeg and meg
NeuroImage 30:1255–1272.

https://doi.org/10.1016/j.neuroimage.2005.10.045
- PubMed
- Google Scholar
1. David O
2. Guillemain I
3. Saillet S
4. Reyt S
5. Deransart C
6. Segebarth C
7. Depaulis A
(2008) Identifying neural drivers with functional MRI: an electrophysiological validation
PLOS Biology 6:e315.

https://doi.org/10.1371/journal.pbio.0060315
- Google Scholar
1. Deshpande G
2. LaConte S
3. James GA
4. Peltier S
5. Hu X
(2009) Multivariate granger causality analysis of fMRI data
Human Brain Mapping 30:1361–1373.

https://doi.org/10.1002/hbm.20606
- PubMed
- Google Scholar
(2008) Analyzing information flow in brain networks with nonparametric granger causality
NeuroImage 41:354–362.

https://doi.org/10.1016/j.neuroimage.2008.02.020
- PubMed
- Google Scholar
(2007) Imaging recollection and familiarity in the medial temporal lobe: a three-component model
Trends in Cognitive Sciences 11:379–386.

https://doi.org/10.1016/j.tics.2007.08.001
- PubMed
- Google Scholar
1. Eden UT
2. Frank LM
3. Barbieri R
4. Solo V
5. Brown EN
(2004) Dynamic analysis of neural encoding by point process adaptive filtering
Neural Computation 16:971–998.

https://doi.org/10.1162/089976604773135069
- PubMed
- Google Scholar
(2007) The medial temporal lobe and recognition memory
Annual Review of Neuroscience 30:123–152.

https://doi.org/10.1146/annurev.neuro.30.051606.094328
- PubMed
- Google Scholar
1. Ekstrom AD
2. Caplan JB
3. Ho E
4. Shattuck K
5. Fried I
6. Kahana MJ
(2005) Human hippocampal theta activity during virtual navigation
Hippocampus 15:881–889.

https://doi.org/10.1002/hipo.20109
- PubMed
- Google Scholar
1. Engle RF
(1982) Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflation
Econometrica: Journal of the Econometric Society 1007:.

https://doi.org/10.2307/1912773
- Google Scholar
1. Engle RF
2. Bollerslev T
(1986) Modelling the persistence of conditional variances
Econometric Reviews 5:1–50.

https://doi.org/10.1080/07474938608800095
- Google Scholar
1. Engle RF
2. Patton AJ
(2001) What good is a volatility model?
Quantitative Finance 1:237–245.

https://doi.org/10.1088/1469-7688/1/2/305
- Google Scholar
1. Ezzyat Y
2. Kragel JE
3. Burke JF
4. Levy DF
5. Lyalenko A
6. Wanda P
7. O'Sullivan L
8. Hurley KB
9. Busygin S
10. Pedisich I
11. Sperling MR
12. Worrell GA
13. Kucewicz MT
14. Davis KA
15. Lucas TH
16. Inman CS
17. Lega BC
18. Jobst BC
19. Sheth SA
20. Zaghloul K
21. Jutras MJ
22. Stein JM
23. Das SR
24. Gorniak R
25. Rizzuto DS
26. Kahana MJ
(2017) Direct brain stimulation modulates encoding states and memory performance in humans
Current Biology 27:1251–1258.

https://doi.org/10.1016/j.cub.2017.03.028
- PubMed
- Google Scholar
1. Ezzyat Y
2. Wanda PA
3. Levy DF
4. Kadel A
5. Aka A
6. Pedisich I
7. Sperling MR
8. Sharan AD
9. Lega BC
10. Burks A
11. Gross RE
12. Inman CS
13. Jobst BC
14. Gorenstein MA
15. Davis KA
16. Worrell GA
17. Kucewicz MT
18. Stein JM
19. Gorniak R
20. Das SR
21. Rizzuto DS
22. Kahana MJ
(2018) Closed-loop stimulation of temporal cortex rescues functional networks and improves memory
Nature Communications 9:365.

https://doi.org/10.1038/s41467-017-02753-0
- PubMed
- Google Scholar
1. Fell J
2. Klaver P
3. Lehnertz K
4. Grunwald T
5. Schaller C
6. Elger CE
7. Fernández G
(2001) Human memory formation is accompanied by rhinal-hippocampal coupling and decoupling
Nature Neuroscience 4:1259–1264.

https://doi.org/10.1038/nn759
- PubMed
- Google Scholar
1. Fell J
2. Fernández G
3. Klaver P
4. Axmacher N
5. Mormann F
6. Haupt S
7. Elger CE
(2006) Rhinal-hippocampal coupling during declarative memory formation: dependence on item characteristics
Neuroscience Letters 407:37–41.

https://doi.org/10.1016/j.neulet.2006.07.074
- PubMed
- Google Scholar
1. Galka A
2. Wong K
3. Ozaki T
(2010)
Modeling Phase Transitions in the Brain

27–52, Generalized state-space models for modeling nonstationary eeg time-series, Modeling Phase Transitions in the Brain, Springer, 10.1007/978-1-4419-0796-7_2.
- Google Scholar
(2014) Understanding predictive information criteria for bayesian models
Statistics and Computing 24:997–1016.

https://doi.org/10.1007/s11222-013-9416-2
- Google Scholar
(2003) Investigating directed cortical interactions in time-resolved fMRI data using vector autoregressive modeling and Granger causality mapping
Magnetic Resonance Imaging 21:1251–1261.

https://doi.org/10.1016/j.mri.2003.08.026
- PubMed
- Google Scholar
1. Goldman RI
2. Stern JM
3. Engel J
4. Cohen MS
(2002) Simultaneous EEG and fMRI of the alpha rhythm
NeuroReport 13:2487–2492.

https://doi.org/10.1097/00001756-200212200-00022
- PubMed
- Google Scholar
Book
(2016)
Deep Learning

MIT Press.
- Google Scholar
1. Graham S
2. Phua E
3. Soon CS
4. Oh T
5. Au C
6. Shuter B
7. Wang S-C
8. Yeh IB
(2009) Role of medial cortical, hippocampal and striatal interactions during cognitive set-shifting
NeuroImage 45:1359–1367.

https://doi.org/10.1016/j.neuroimage.2008.12.040
- Google Scholar
(1994) Multivariate stochastic variance models
The Review of Economic Studies 61:247–264.

https://doi.org/10.2307/2297980
- Google Scholar
1. Harvey AC
2. Shephard N
(1996) Estimation of an asymmetric stochastic volatility model for asset returns
Journal of Business & Economic Statistics 14:429–434.

https://doi.org/10.1080/07350015.1996.10524672
- Google Scholar
1. Hasselmo ME
(2005) What is the function of hippocampal theta rhythm?--linking behavioral data to phasic properties of field potential and unit recording data
Hippocampus 15:936–949.

https://doi.org/10.1002/hipo.20116
- PubMed
- Google Scholar
(2004) Cognitive functions of gamma-band activity: memory match and utilization
Trends in Cognitive Sciences 8:347–355.

https://doi.org/10.1016/j.tics.2004.06.006
- PubMed
- Google Scholar
1. Hesse W
2. Möller E
3. Arnold M
4. Schack B
(2003) The use of time-variant EEG Granger causality for inspecting directed interdependencies of neural assemblies
Journal of Neuroscience Methods 124:27–44.

https://doi.org/10.1016/S0165-0270(02)00366-7
- PubMed
- Google Scholar
1. Heston SL
(1993) A Closed-Form solution for options with stochastic volatility with applications to bond and currency options
Review of Financial Studies 6:327–343.

https://doi.org/10.1093/rfs/6.2.327
- Google Scholar
1. Heston SL
2. Nandi S
(2000) A Closed-Form GARCH option valuation model
Review of Financial Studies 13:585–625.

https://doi.org/10.1093/rfs/13.3.585
- Google Scholar
Book
1. Hull J
(2009)
Options, Futures and Other Derivatives

Pearson/Prentice Hall.
- Google Scholar
1. Jacobs J
2. Kahana MJ
(2010) Direct brain recordings fuel advances in cognitive electrophysiology
Trends in Cognitive Sciences 14:162–171.

https://doi.org/10.1016/j.tics.2010.01.005
- PubMed
- Google Scholar
(1994) Bayesian analysis of stochastic volatility models
Journal of Business & Economic Statistics 12:371–389.

https://doi.org/10.1080/07350015.1994.10524553
- Google Scholar
1. Jensen O
2. Colgin LL
(2007) Cross-frequency coupling between neuronal oscillations
Trends in Cognitive Sciences 11:267–269.

https://doi.org/10.1016/j.tics.2007.05.003
- PubMed
- Google Scholar
(2001) Theta returns
Current Opinion in Neurobiology 11:739–744.

https://doi.org/10.1016/S0959-4388(01)00278-1
- PubMed
- Google Scholar
1. Kass RE
2. Raftery AE
(1995) Bayes factors
Journal of the American Statistical Association 90:773–795.

https://doi.org/10.1080/01621459.1995.10476572
- Google Scholar
1. Kastner G
2. Frühwirth-Schnatter S
(2014) Ancillarity-sufficiency interweaving strategy (ASIS) for boosting MCMC estimation of stochastic volatility models
Computational Statistics & Data Analysis 76:408–423.

https://doi.org/10.1016/j.csda.2013.01.002
- Google Scholar
1. Kelley WM
2. Miezin FM
3. McDermott KB
4. Buckner RL
5. Raichle ME
6. Cohen NJ
7. Ollinger JM
8. Akbudak E
9. Conturo TE
10. Snyder AZ
11. Petersen SE
(1998) Hemispheric specialization in human dorsal frontal cortex and medial temporal lobe for verbal and nonverbal memory encoding
Neuron 20:927–936.

https://doi.org/10.1016/S0896-6273(00)80474-2
- PubMed
- Google Scholar
1. Kim S
2. Shepherd N
3. Chib S
(1998) Stochastic volatility: likelihood inference and comparison with ARCH models
Review of Economic Studies 65:361–393.

https://doi.org/10.1111/1467-937X.00050
- Google Scholar
1. Kirwan CB
2. Stark CE
(2004) Medial temporal lobe activation during encoding and retrieval of novel face-name pairs
Hippocampus 14:919–930.

https://doi.org/10.1002/hipo.20014
- PubMed
- Google Scholar
1. Klimesch W
(1999) EEG alpha and theta oscillations reflect cognitive and memory performance: a review and analysis
Brain Research Reviews 29:169–195.

https://doi.org/10.1016/S0165-0173(98)00056-3
- PubMed
- Google Scholar
(2004) ‘Epileptic seizures’ in economic organism
Physica A: Statistical Mechanics and Its Applications 342:701–707.

https://doi.org/10.1016/j.physa.2004.05.045
- Google Scholar
1. Kreiman G
2. Koch C
3. Fried I
(2000) Category-specific visual responses of single neurons in the human medial temporal lobe
Nature Neuroscience 3:946–953.

https://doi.org/10.1038/78868
- PubMed
- Google Scholar
(2014) Cross-validation pitfalls when selecting and assessing regression and classification models
Journal of Cheminformatics 6:10.

https://doi.org/10.1186/1758-2946-6-10
- PubMed
- Google Scholar
(1996) Effect of age on the sleep EEG: slow-wave activity and spindle frequency activity in young and middle-aged men
Brain Research 738:205–212.

https://doi.org/10.1016/S0006-8993(96)00770-6
- PubMed
- Google Scholar
(2003) EEG-correlated fMRI of human alpha activity
NeuroImage 19:1463–1476.

https://doi.org/10.1016/S1053-8119(03)00286-6
- PubMed
- Google Scholar
(2014) Subsequent memory effect in intracranial and scalp EEG
NeuroImage 84:488–494.

https://doi.org/10.1016/j.neuroimage.2013.08.052
- PubMed
- Google Scholar
1. Long NM
2. Kahana MJ
(2015) Successful memory formation is driven by contextual encoding in the core memory network
NeuroImage 119:332–337.

https://doi.org/10.1016/j.neuroimage.2015.06.073
- PubMed
- Google Scholar
1. Luo Q
2. Ge T
3. Grabenhorst F
4. Feng J
5. Rolls ET
(2013) Attention-dependent modulation of cortical taste circuits revealed by Granger causality with signal-dependent noise
PLoS Computational Biology 9:e1003265.

https://doi.org/10.1371/journal.pcbi.1003265
- PubMed
- Google Scholar
1. Macke JH
2. Buesing L
3. Cunningham JP
4. Byron MY
5. Shenoy KV
6. Sahani M
(2011)
Empirical models of spiking in neural populations

Advances in Neural Information Processing Systems.
- Google Scholar
1. Maloney KJ
2. Cape EG
3. Gotman J
4. Jones BE
(1997) High-frequency gamma electroencephalogram activity in association with sleep-wake states and spontaneous behaviors in the rat
Neuroscience 76:541–555.

https://doi.org/10.1016/S0306-4522(96)00298-9
- PubMed
- Google Scholar
1. Manning JR
2. Jacobs J
3. Fried I
4. Kahana MJ
(2009) Broadband shifts in local field potential power spectra are correlated with single-neuron spiking in humans
Journal of Neuroscience 29:13613–13620.

https://doi.org/10.1523/JNEUROSCI.2041-09.2009
- PubMed
- Google Scholar
(2017) Arima-garch 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
- Google Scholar
1. Nelson DB
(1991) Conditional heteroskedasticity in asset returns: a new approach
Econometrica: Journal of the Econometric Society.

https://doi.org/10.2307/2938260
- Google Scholar
1. Nir Y
2. Staba RJ
3. Andrillon T
4. Vyazovskiy VV
5. Cirelli C
6. Fried I
7. Tononi G
(2011) Regional slow waves and spindles in human sleep
Neuron 70:153–169.

https://doi.org/10.1016/j.neuron.2011.02.043
- PubMed
- Google Scholar
1. Ojemann GA
2. Dodrill CB
(1985) Verbal memory deficits after left temporal lobectomy for epilepsy. mechanism and intraoperative prediction
Journal of Neurosurgery 62:101–107.

https://doi.org/10.3171/jns.1985.62.1.0101
- PubMed
- Google Scholar
1. Omori Y
2. Chib S
3. Shephard N
4. Nakajima J
(2007) Stochastic volatility with leverage: fast and efficient likelihood inference
Journal of Econometrics 140:425–449.

https://doi.org/10.1016/j.jeconom.2006.07.008
- Google Scholar
1. Pedregosa F
2. Varoquaux G
3. Gramfort A
4. Michel V
5. Thirion B
6. Grisel O
7. Blondel M
8. Prettenhofer P
9. Weiss R
10. Dubourg V
11. Vanderplas J
12. Passos A
13. Cournapeau D
14. Brucher M
15. Perrot M
16. Duchesnay E
(2011)
Scikit-learn: machine learning in Python

Journal of Machine Learning Research 12:2825–2830.
- Google Scholar
(2001) Gating of human theta oscillations by a working memory task
The Journal of Neuroscience 21:3175–3183.

https://doi.org/10.1523/JNEUROSCI.21-09-03175.2001
- PubMed
- Google Scholar
1. Rasmussen CE
(2004)
Advanced Lectures on Machine Learning

63–71, Gaussian processes in machine learning, Advanced Lectures on Machine Learning, Springer, 10.1007/978-3-540-28650-9_4.
- Google Scholar
(2005) Mapping directed influence over the brain using Granger causality and fMRI
NeuroImage 25:230–242.

https://doi.org/10.1016/j.neuroimage.2004.11.017
- PubMed
- Google Scholar
(2003) Theta and gamma oscillations during encoding predict subsequent recall
The Journal of Neuroscience 23:10809–10814.

https://doi.org/10.1523/JNEUROSCI.23-34-10809.2003
- PubMed
- Google Scholar
(2007) Hippocampal and neocortical gamma oscillations predict memory formation in humans
Cerebral Cortex 17:1190–1196.

https://doi.org/10.1093/cercor/bhl030
- PubMed
- Google Scholar
(2015) Granger causality analysis in neuroscience and neuroimaging
Journal of Neuroscience 35:3293–3297.

https://doi.org/10.1523/JNEUROSCI.4399-14.2015
- PubMed
- Google Scholar
1. Smith AC
2. Brown EN
(2003) Estimating a State-Space model from point process observations
Neural Computation 15:965–991.

https://doi.org/10.1162/089976603765202622
- Google Scholar
1. Solomon EA
2. Kragel JE
3. Sperling MR
4. Sharan A
5. Worrell G
6. Kucewicz M
7. Inman CS
8. Lega B
9. Davis KA
10. Stein JM
11. Jobst BC
12. Zaghloul KA
13. Sheth SA
14. Rizzuto DS
15. Kahana MJ
(2017) Widespread theta synchrony and high-frequency desynchronization underlies enhanced cognition
Nature Communications 8:1704.

https://doi.org/10.1038/s41467-017-01763-2
- PubMed
- Google Scholar
1. Solomon EA
2. Stein JM
3. Das S
4. Gorniak R
5. Sperling MR
6. Worrell G
7. Inman CS
8. Tan RJ
9. Jobst BC
10. Rizzuto DS
11. Kahana MJ
(2019) Dynamic theta networks in the human medial temporal lobe support episodic memory
Current Biology 29:1100–1111.

https://doi.org/10.1016/j.cub.2019.02.020
- PubMed
- Google Scholar
(2002) Bayesian measures of model complexity and fit
Journal of the Royal Statistical Society: Series B 64:583–639.

https://doi.org/10.1111/1467-9868.00353
- Google Scholar
1. Squire LR
2. Zola-Morgan S
(1991) The medial temporal lobe memory system
Science 253:1380–1386.

https://doi.org/10.1126/science.1896849
- PubMed
- Google Scholar
(2010) Ten simple rules for dynamic causal modeling
NeuroImage 49:3099–3109.

https://doi.org/10.1016/j.neuroimage.2009.11.015
- PubMed
- Google Scholar
1. Stokes PA
2. Purdon PL
(2017) A study of problems encountered in Granger causality analysis from a neuroscience perspective
PNAS 114:E7063–E7072.

https://doi.org/10.1073/pnas.1704663114
- PubMed
- Google Scholar
1. Stone M
(1977) An asymptotic equivalence of choice of model by cross-validation and akaike’s criterion
Journal of the Royal Statistical Society. Series B 39:44–47.

https://doi.org/10.1111/j.2517-6161.1977.tb01603.x
- Google Scholar
1. von Stein A
2. Sarnthein J
(2000) Different frequencies for different scales of cortical integration: from local gamma to long range alpha/theta synchronization
International Journal of Psychophysiology 38:301–313.

https://doi.org/10.1016/S0167-8760(00)00172-0
- PubMed
- Google Scholar
1. Wang Y
(2002) Asymptotic nonequivalence of GARCH models and diffusions
The Annals of Statistics 30:754–783.

https://doi.org/10.1214/aos/1028674841
- Google Scholar
Book
1. West M
(1996)
Bayesian Forecasting

Wiley Online Library.
- Google Scholar
1. Wong KF
2. Galka A
3. Yamashita O
4. Ozaki T
(2006) Modelling non-stationary variance in EEG time series by state space GARCH model
Computers in Biology and Medicine 36:1327–1335.

https://doi.org/10.1016/j.compbiomed.2005.10.001
- PubMed
- Google Scholar
1. Wu W
2. Gao Y
3. Bienenstock E
4. Donoghue JP
5. Black MJ
(2006) Bayesian population decoding of motor cortical activity using a kalman filter
Neural Computation 18:80–118.

https://doi.org/10.1162/089976606774841585
- PubMed
- Google Scholar
1. Yu Y
2. Meng X-L
(2011) To center or not to center: that is not the Question—An Ancillarity–Sufficiency Interweaving Strategy (ASIS) for Boosting MCMC Efficiency
Journal of Computational and Graphical Statistics 20:531–570.

https://doi.org/10.1198/jcgs.2011.203main
- Google Scholar
1. Yushkevich PA
2. Pluta JB
3. Wang H
4. Xie L
5. Ding SL
6. Gertje EC
7. Mancuso L
8. Kliot D
9. Das SR
10. Wolk DA
(2015) Automated volumetry and regional thickness analysis of hippocampal subfields and medial temporal cortical structures in mild cognitive impairment
Human Brain Mapping 36:258–287.

https://doi.org/10.1002/hbm.22627
- PubMed
- Google Scholar
1. Zhou Z
2. Ding M
3. Chen Y
4. Wright P
5. Lu Z
6. Liu Y
(2009) Detecting directional influence in fMRI connectivity analysis using PCA based Granger causality
Brain Research 1289:22–29.

https://doi.org/10.1016/j.brainres.2009.06.096
- PubMed
- Google Scholar

Article and author information

Author details

Tung D Phan

University of Pennsylvania, Philadelphia, United States

Present address
Department of Psychology, University of Pennsylvania, Philadelphia, United States

Contribution
Conceptualization, Data curation, Software, Formal analysis, Validation, Visualization, Writing—original draft, Writing—review and editing

For correspondence
tungphan87@gmail.com

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-5957-7566
Jessica A Wachter

University of Pennsylvania, Philadelphia, United States

Present address
Department of Finance, The Wharton School, University of Pennsylvania, Philadelphia, United States

Contribution
Conceptualization, Supervision, Validation, Writing—review and editing

Competing interests
No competing interests declared
Ethan A Solomon

University of Pennsylvania, Philadelphia, United States

Present address
Department of Bioengineering, University of Pennsylvania, Philadelphia, United States

Contribution
Writing—review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-0541-7588
Michael J Kahana

University of Pennsylvania, Philadelphia, United States

Present address
Department of Psychology, University of Pennsylvania, Philadelphia, United States

Contribution
Conceptualization, Supervision, Funding acquisition, Investigation, Methodology, Writing—original draft, Writing—review and editing

For correspondence
kahana@psych.upenn.edu

Competing interests
No competing interests declared

Funding

Defense Advanced Research Projects Agency (N66001-14-2-4032)

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 N66001-14-2-4032). 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, Dartmouth-Hitchcock 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.

Copyright

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.