Extracting single-trial neural interaction using latent dynamical systems model

In systems neuroscience, advances in simultaneous recording technology have helped reveal the population dynamics that underlie the complex neural correlates of animal behavior and cognitive processes. To investigate these correlates, neural interactions are typically abstracted from spike trains of pairs of neurons accumulated over the course of many trials. However, the resultant averaged values do not lead to understanding of neural computation in which the responses of populations are highly variable even under identical external conditions. Accordingly, neural interactions within the population also show strong fluctuations. In the present study, we introduce an analysis method reflecting the temporal variation of neural interactions, in which cross-correlograms on rate estimates are applied via a latent dynamical systems model. Using this method, we were able to predict time-varying neural interactions within a single trial. In addition, the pairwise connections estimated in our analysis increased along behavioral epochs among neurons categorized within similar functional groups. Thus, our analysis method revealed that neurons in the same groups communicate more as the population gets involved in the assigned task. We also showed that the characteristics of neural interaction from our model differ from the results of a typical model employing cross-correlation coefficients. This suggests that our model can extract nonoverlapping information about network topology, unlike the typical model.


Introduction
Information communication via spike trains of neurons in populations is a core computational process that enables many brain areas to execute their roles, which include encoding of stimuli, decision-making, and highlevel cognition [1]. To understand these processes, therefore, the effects of spike trains across neuronal populations must be determined according to their specific network structures [2][3][4][5][6]. For this purpose, a variety of network theoretical models have been developed for the analysis of neuronal population dynamics and the description of network topology, e.g., modules, hubs, and rich-clubs [7][8][9]. From analyses of pairwise correlations, a recent study showed that the functional single neuron network of the macaque monkey frontoparietal area during active behavior has a highly complex topology that includes small-worldness, hubs, and rich-clubs consisting of oscillatory neurons synchronized in specific frequency bands [9]. To infer neural network topology, many studies have focused on quantifying correlations in the number of spikes or in the spike times of pairwise neurons. One conventionally used quantification method is the cross-correlogram (CCG) coefficient, which measures the co-occurrence of spikes of the pairs within a given time bin [10]. To assess temporal correlations rather than the degree of synchronous firing, alternative quantification methods can be applied, such as the correlation index [11] or spike-time tiling coefficient [12], which are derived from the temporal order of pairwise spike trains. Given sufficient data, it is also possible to make use of nonlinear methods such as mutual information [13] and transfer entropy [14] in information theory.
Measuring sets of neural responses across trials is common to estimate the firing rates of spike trains because spikes are highly variable and noisy at the single-trial level. The averaged values across trials, however, obfuscate information of population dynamics at certain times in behavioral epochs or in the trials. Since estimating neural activities in single trials is challenging because of variability and deficient sampling, several studies have exploited the dimensionality reduction method adapted in the latent variable space rather than the neuronal space to measure population activities at the single-trial level [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. Models with independent factors underlying population activities can extract neural population states in single trials [15,16]. Recent studies have reported analysis of single-trial population activities using linear dynamical models across time [17][18][19][20] or nonlinear dynamical models allowing switches of states at given times [21,22]. Event-dependent dynamical systems models [22] are variant models compromising linear and nonlinear methods that allow switches at the cued times between adjacent epochs. In addition, artificial recurrent neural networks infer single-trial neural population dynamics based on the assumption that the networks can generate neural data using a machinelearning method [23]. In nonhuman primates, population activities in diverse brain regions are low dimensional in many instances [19,[30][31][32][33][34][35][36][37]; therefore, it can be assumed that population activity in single trials can be estimated via a latent dynamical systems model with much lower dimensions than the population size.
Estimations of population activity are useful for elucidating the slow-varying neural dynamics governed by shared inputs in a single trial; however, the recovery of transient changes evoked by other neurons' spikes is limited because the temporal precision of the latent model is generally insufficient for capturing these interactions. The intrinsic timescale of spike effects is about tens of ms and below 200 ms at most [38], whereas neural populations sustain their activity in the same-shared activity space within an epoch that generally lasts longer than hundreds of ms [22]. To address this issue, it is possible to calculate CCGs of pairs of neurons by two different ways that CCGs are calculated using spike trains or rate estimates. In typical fashion, CCGs are calculated using spike trains, but they include the effects caused by shared motives such as common inputs or synchronization across the populations as well as the direct effects of the other neurons in the population. To quantify these shared terms, the CCGs can be calculated using estimates of population activity inferred by a latent dynamical systems model. The transient effects between pairs of neurons can be surmised by comparing two CCGs according to spike trains or estimates of populate activity. Estimates of population activity can then be revised in light of transient effects, and the degree to which both transient and shared effects contribute to the revised estimate can be evaluated.
In the present study, we applied the above-described method to a dataset from an open database in which data were simultaneously recorded from the anterior lateral motor (ALM) cortex of mice while they were executing a two-alternative-forced choice task [29,[39][40][41][42][43]. We found that neural interactions varied significantly over time in a trial and that neurons could be classified based on the degree of other neurons' effects in each behavioral epoch. We also showed that network structures were dependent on the task-related epochs, implying that interactions between neurons were not stationary but instead adjusted to the relative importance of the epochs. Compared with a model containing correlation coefficients, which are widely used to evaluate functional connectivity, our model could improve judgment of the underlying organization of transient neural activities.

Application to mice ALM cortex neurons
We applied the above-described analysis method to a dataset from alternative choice tasks conducted in previous studies [29,[39][40][41][42][43]. In this research, trained mice reported pole position (posterior or anterior) by licking one of two targets (left in an ipsi trial and right in a contra trial) following a delay epoch after pole presentation. After a 1.3-s delay, mice responded to an auditory signal indicating the onset of the response epoch (Fig. 1a, upper panel). Extracellular spikes were recorded on the left-hemisphere ALM cortex using 32-channel Neu-roNexus silicon probes (19 mice) or 64-channel Janelia silicon probes (20 mice). The extracellular recording traces were bandpass-filtered (300-6000 Hz) and then sorted using JRclust if events exceeded an amplitude threshold. Based on spike width, cells were classified as pyramidal cells (width > 0.5 ms) or fast-spiking interneurons (width < 0.35 ms). For a more detailed explanation of the electrophysiological recording method, consult the studies cited earlier [40,42,43]. To assess the effects of ensemble inputs from other neurons, we analyzed only the sessions that included > 15 neurons (9 of 40 sessions).

Estimating the dynamics of ALM cortex neural activities
Neural dynamics are highly complicated and stochastic, resulting in diverse responses even in comparable cases. To resolve these difficulties, the extraction of shared neural activity (SNA) in a low-dimensional latent space has been developed to determine the dynamics of population activity as an alternative to directly assessing raw neural activities. For this purpose, we used a previously introduced model [22], i.e., an event-dependent linear dynamical systems (EDLDS) model that allows switches of behavioral epochs as follows: The latent variable x(t) is an N L -dimensional vector ( N L : the dimension of the latent space) updated by matrix W (s) in the s-epoch (s: sample, delay, response). The term W 0 (t) represents random fluctuation independent of the term updated by matrix W (s) in the latent space. The rate estimate r(t) is an N P -dimensional vector ( N P : the number of neurons in the population) calculated by the projection matrix V (s) in the s-epoch. The (1) term V 0 (t) denotes the residual that is not explained by the latent variables.
Temporal interaction estimated with CCGs. Since the EDLDS model captures slow-varying dynamics of the neural response, neural components that reflect temporal fluctuations must be introduced into the model. One such possible component is the interaction between neurons over a finite time scale. We assumed that a neuron was affected by all the other neurons in the population and that the rate estimate calculated in Eq. 1 could be revised by the integration of these interactions. Accordingly, we calculated CCGs between all pairs of neurons from their spike trains C ij (kτ ) as follows: where j ( i ) is the index of a reference (a target) neuron affecting (affected by) the other neuron. M and N m j are n m i t m jl , k /τ of the reference neuron. We found similar results with various temporal precisions from τ = 15ms to τ = 45ms for the transient effects; thus, we heuristically selected the temporal precision τ = 25ms because the results of analysis using this value showed greater distinction in neural interactions according to the epoch. In a similar manner, the estimated CCG (ECCG) from the SNA of a latent model C E ij (kτ ) can be represented as: where r m i t m jl , k denotes the mean rate estimate of the target neuron in the k-th time lag range t m jl + (k − 1)τ ≤ t < t m jl + kτ after the l-th spike time t m jl of the reference neuron. In line with n m i t m jl , k /τ in Eq. 3, which denotes the rate in the range t m jl + (k − 1)τ ≤ t < t m jl + kτ derived from the spike trains of the i-th neuron, the counterpart r m i t m jl , k in the ECCG (Eq. 4) is the rate in the same range derived from the rate estimate of the i-th neuron (Eq. 1). The difference between the CCG and ECCG (i.e., the DCCG) �C ij (kτ ) can then be represented as follows:

Optimization of neural activities
SNA effectively describes the slow-varying population activity governed by the shared inputs. However, recovering the transient changes evoked by the spikes of the connected neurons is difficult due to the short time variability of the interaction effect compared with the temporal precision of the latent model. For increasingly accurate estimates that include these transient effects over a relatively short time scale, the rate estimate in the given time bin of the latent model can be corrected using the CCG and ECCG. Since the interaction between neurons might change across times or trials, we allowed the CCG and ECCG to depend on both variables. We calculated the CCG and ECCG from the spike trains in the neighboring 20 bins of 20 trials (we obtained similar results using different values, e.g., 15 or 30 bins or trials; data not shown). For simplicity, we used only one-time lag before and after the spike of the reference neuron ( C ij (τ ), C ij (−τ ) ) and then the corrected rate estimate r T i,C (t) as follows: , T denotes the time bin index of the latent model, and 0 ≤ t < S bin ( S bin : time bin size of the latent model). We used S bin = 67.4 ms following a previous study [22]. The rate estimate r T i represents the SNA of the i-th neuron in the given time bin T . Because of the effects of the other neurons' spikes, the mean of the corrected rate estimate r T i,C (t) in the given time bin T becomes different from the rate estimate r T i . To resolve this disagreement, we normalized the corrected rate estimate to match with the SNA; the transient SNA (TSNA) r T i,TSNA (t) was then represented as follows: where r T i,C (t) is the mean of r T i,C (t) in the given time bin T . Since the contributions of transient effects to the neural response vary depending on the time bin, we found the optimal degree of SNA and TSNA contributions by maximizing the likelihood L({t i }) of the target neuron's spike times in a given time bin as follows: where {t i } denotes the spike trains of neuron i and δ t,{t i } = 1 or 0 if t ∈ {t i } or does not, respectively. The probability P(t) of firing at time t is calculated using the soft-max selection rule . We set a discrete time step of δt = S bin /60(S bin = 67.4 ms: time bin size of the latent model) for all calculations in this work. The parameter α denotes the degree to which TSNA contributes and it is determined by Q. The optimized neural activity (ONA)r T i,ONA (t) is then represented as follows: where α * denotes the value maximizing L({t i }) , i.e., α * = 1/(1 + e −Q * ) and Q * = argmax Q L({t i }).

Neural interactions estimated by Kullback-Leibler divergence
To assess how SNA varies by the transient effects of other neurons, we compared the probability distributions of the spike trains, i.e., Kullback-Leibler divergences, between the ONA and SNA (i.e., the DONA) D KL T (P ONA T (t)|P SNA T ) as follows: where P SNA T ( P ONA T (t) ) denotes the probability of firing for SNA(ONA) at time t in the given time bin T . To identify interactions induced by a single reference neuron's spike train but not by all other neurons (as with DONA), we defined the pairwise Kullback-Leibler divergence between ONA and SNA (i.e., the PDONA) D KL T P ONA (s n )|P SNA (s n ) as follows: where s n is an n-dimensional vector with 0 or 1 denoting spike existence or nonexistence of the target neuron i in the time step δt after n-consecutive spike times of the reference neuron j ( S(n) : the entire space including all the possible cases of s n ) closest to the time bin T, whereas t k j denotes the k-th spike time of the reference neuron j, and δ 1,s n (k) = 1 or 0 if the k-th element of s n ( s n (k) ) is 1 or 0, respectively. P(s n ) represents the probability of the series of the target neuron firing or not firing after the firing of the reference neuron. Hence, PDONA measures how each of the reference neuron's n-consecutive spikes modify the probability of the target neuron firing after the moment at which the corresponding reference neuron spikes. We tested various cases of n spikes of reference neurons (from n = 5 to n = 8); we chose to use n = 6 and excluded cases in which the number of spikes by the reference neuron was < 10 in the trial (we obtained similar results for n = 5, 7, or 8; data not shown).

Logistic regression analysis for decoding lick responses
We evaluated how the neural interaction improved performance in decoding lick responses through logistic regression analysis, including spike counts and DONA (full model), as follows: where the log-odds for the probability of the lick of the left target in the ipsi trial ( P ipsi ) were fitted using the log of the spike counts ( log(N sp ) ) and the log of DONA ( log(DONA) ). By comparing Akaike's information criterion (AIC) between the full model and reduced model including only the spike counts, we determined the optimal model for decoding lick responses.

Cross-correlation coefficients
To compare DONA with a typical measure of neural interaction, we calculated the cross-correlation coefficients of neuron pairs in the given time bin. For direct comparison with DONA, we defined a normalized correlation coefficient (NCC) ρ i , which assumes the mean correlation between a single neuron i and all the other neurons as follows: (13) log P ipsi /(1 − P ipsi ) where N P is the number of neurons in the population and ρ ij is the cross-correlation coefficient between neuron i and neuron j.

Population activities and CCGs
To demonstrate the neural interactions in an ALM cortex, we analyzed neural data collected from mice as they executed a delayed response task (Fig. 1a, upper panel). ALM neurons showed complex and variable activities, as shown in the example raster plot and firing rates from a single trial (Fig. 1a, middle and   and ECCGs were sustained over hundreds of ms as they moved farther away from the reference neuron's spike time (time = 0). DCCGs, however, approached zero after a few bins (about tens of ms), implying that the DCCG captures the transient effect near time zero by offsetting the shared effects of both the CCG and ECCG remaining in the longer time range.

Localized neural interaction terms at specific behavioral epochs
In our analysis, DONA is a quantity with which we can judge how the probability distribution of the ONA differs from that of the SNA, which is acquired from the latent model ( Eqs. 1-2). Figure 2a provides an example of the optimization of the activities of the target neuron N15 (Fig. 1b), including all the effects of the other (reference) neurons, with parameter α deciding the degree of SNA and TSNA contributions (Fig. 2b). Since ONA is the optimal neural activity between the transient component from other neurons and shared components due to behavior, larger values of DONA can be interpreted as more interactions occurring with the other neurons. Neuron N15 in Fig. 2c, for example, shows more interactions with other neurons during the sensory (indicated by S) and response (indicated by R) epochs compared with the interaction observed during the delay epoch (indicated by D). For pairwise interactions, we showed PDO-NAs of three pairs, which represent the effects of three reference neurons N2, N7, and N12 on the target neuron N15 shown in Fig. 1b (Fig. 2d). It should be noted that the DONA of the target neuron N15 in Fig. 2c was the interaction integrated with all other neurons, although we only showed the results of the three reference neurons to improve visibility. To characterize the neural interactions of all neurons, we classified them in conformity with the epochs in which their DONAs were maximized. First,  Fig. 3a (session no. 37), DONAs were significantly larger; therefore, they were well localized in specific epochs. Among 22 neurons, 6 (27%), 5 (23%), and 5 (23%) showed more interactions with other neurons in the sample, delay and response behavioral epochs, respectively. For the neurons of all sessions (157 neurons in 9 sessions), most (111 of 157; 70.70%) showed maximum values of DONA consistently in single epochs (Fig. 3b). This result implies that the effects of neural interactions on single neurons are rather concentrated on a moment of information processing than on being maintained throughout the entire trial. The remaining neurons showed larger values of DONA during the two epochs (26 of 157; 16.54%) or no larger values (20 of 157; 12.74%).

Tendency of increases of neural interactions toward the end of a trial
To determine the overall variation of DONA during the progress of epochs, we simply averaged DONA values for all the neurons regardless of their statistical significance in epochs (Fig. 3c). The mean DONA value increased and reached a peak soon after the response onset. Significantly more neurons showed larger DONA values in the latter epochs, in agreement with the results of mean DONA values (p < 0.05, χ 2 test; Fig. 3d). The number of neurons with larger DONA values was 33, 50, and 68 of 157 total neurons during the sample, delay, and response epochs, respectively; thus, there were significant differences between the behavioral epochs (p < 0.05 between S and D and between D and R, both χ 2 tests).
In the dataset used for this analysis, Hidehiko et al. [43] categorized neurons in the ALM cortex into four functional groups, namely contralateral ramping-up, ipsilateral ramping-up, ramping-down, and the rest based on the selectivity of the lick responses and the neurons' ramping activities. We investigated whether the neuronal interaction term DONA was related to these functional neuron groups. According to the categorization of Hidehiko et al. [43], we placed 52, 17, 19, and 69 neurons into the contralateral ramping-up, ipsilateral rampingup, ramping-down, and remaining groups, respectively. Secondly, we calculated PDONAs between neurons belonging to the same functional groups and the different functional groups separately. As was the case for DONA, pairs of neurons were considered to show larger interactions in an epoch if the magnitudes of PDONA in the epoch (20 bins per epoch) were larger than the magnitudes in the other two epochs (Wilcoxon rank sum test, p < 0.05). The proportion of pairs with larger values of PDONA in each epoch is shown Fig. 4a (filled and empty bars show pairs of neurons in the same functional group and different groups, respectively). Significantly more pairs in the same functional group showed larger PDONA values in the latter epochs in accordance with the DONA results shown in Fig. 3d (Fig. 4a, black bars). In the same functional group, the pairs of neurons showing larger PDONA values were 36, 73, and 121 of 184 total pairs during the sample, delay, and response epochs, respectively (p < 0.01 between S and D and between D and R, both χ 2 tests). However, for pairs of neurons belonging to different functional groups, there was no significant difference in the number of pairs with larger PDONA values between epochs (Fig. 4a, empty bars): 374, 403, and 437 of 1,125 total pairs during the sample, delay, and response epochs, respectively (p = 0.199 between S and D; p = 0.138 between D and R; both χ 2 tests). These results imply that the increase in neural interactions in the latter epoch was due to increased functional connections between neurons in the same functional group.
The number of pairs of larger PDONA values was significantly different in the sample and response epochs but not in the delay epoch (delay: p = 0.314, χ 2 test). In the sample epoch, the proportion of pairs with larger PDONA values was significantly higher when neurons in the pairs belonged to different functional groups (p < 0.01, χ 2 test). However, the opposite was true in the response epoch in which the proportion of pairs of neurons in the same functional group was higher (p < 0.01, χ 2 test). We also reanalyzed the data for each group separately (Fig. 4b). The numbers of pairs of neurons in the contralateral ramping-up, ipsilateral ramping-up, and ramping-down group were 150, 18, and 16 of 184 total pairs of neurons. It should be noted that we did not analyze the fourth group (the rest) separately because neurons in the fourth group were just the remainders that could not be classified into any of the three groups discussed above. Although the results for the contralateral ramping-up group were identical to those of all the groups shown in Fig. 4a (Fig. 4b, left panel), the PDONAs of the other groups had different distributions when compared with the overall tendency to increase along the epoch (Fig. 4b, middle and right panels). In the contralateral ramping-up group, the pairs of neurons showing larger PDONA values were 26, 54, 117 of 150 total pairs during the sample, delay, and response epochs, respectively (p < 0.01 between S and D and between D and R, both χ 2 tests). In the ipsilateral ramping-up group, the pairs of neurons showing larger PDONA values were 2, 13, 4 of 18 total pairs during the sample, delay, and response epochs, respectively (p < 0.01 between S and D and between D and R, both χ 2 tests). In the ramping-down group, the pairs of neurons showing larger PDONA values were 8, 6, 0 of 16 total pairs during the sample, delay, and response epochs, respectively (p = 0.48 between S and D; p < 0.01 between D and R; both χ 2 tests). This suggests that neurons classified in different groups show different response patterns. In addition, the interactions between these neurons also show different tendencies according to epoch.

Enhanced decoding capability for behavioral responses due to neural interactions
We performed logistic regression analyses with the spike count N sp and DONA to test whether neural interactions enhanced the performance when decoding lick responses ( Fig. 5a; see Material and methods). The probability of the response was significantly dependent (p < 0.05 in the logistic regression) on DONA in 10.53% of neurons (p < 0.05, binomial test) and on the spike count in 71.58% of neurons (p < 0.01, binomial test), as shown in Fig. 5a. We compared the full model including both DONA and spike counts, which we expected to be optimal for single neurons, with the reduced model including only spike counts. Comparing the AICs of the two models, we found that the full model was optimal in a considerable number of neurons (21.13%; Fig. 5b). Thus, it seems that DONA can convey more information than spike count alone. In Fig. 5c, d, we compare the dynamics of DONA and spike counts according to lick responses by showing a representative example of a single neuron for which the probability of the response depended significantly on both DONA and the spike counts of the neuron in the delay epoch.

Disparate trends of cross-correlation compared with DONA
To compare the localization trends of DONA (shown in Fig. 3a-c) with an existing and typical measure of pairwise interactions, we calculated NCCs for each neuron ( Fig. 6; see Materials and methods). To enable the direct comparison of these measures, the results of NCC were arranged for the neurons sorted in the same manner shown in Fig. 3a-c. If NCC values could extract similar features of network topology as DONA, we would see a resemblance between NCC and DONA values in terms of localization to epochs. However, NCCs did not show such a localization tendency and the mean NCC showed no patterns along the epochs, unlike the mean DONA (Fig. 6b, open circle). These results suggest that DONA might extract nonoverlapping features related to the neural interactions in the population that cannot be extracted by the typically used correlation measure.

Discussion
In the present study, we present a method by which to extract single-trial neural interactions using rate estimates obtained by a latent dynamical model. We applied this method to a dataset containing data from mice cortical neurons [22,[39][40][41][42][43]. The model applied in this method was able to evaluate the proportions of shared neural activities and transient effects evoked by other neurons; thus, it showed that inferred measures could trace the variation of neural interactions along epochs. Most neurons were maximally affected by other neurons in a specific epoch and, on average, neural interactions increased as an epoch progressed. To demonstrate the novelty of our method, we compared its performance to that achieved using cross-correlation coefficients; we showed that our model could reveal hidden features of the neural network that were not accounted for by a typical method in which pairwise interactions are estimated. Obstacles must be overcome to ensure that our method can be generally applied to neural population data in different regions of the brain when performing diverse tasks. First, a number of existing models [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29] have been established by which to analyze single-trial neural population dynamics and provide different rate estimates to certain extents. It might also be more challenging to estimate neural activities in single trials for tasks that require higher computational capabilities [21], although many studies have shown that population activities could be reduced to low-dimensional shared spaces [25,[30][31][32][33][34][35][36][37]. In such tasks, latent dynamical models provide suboptimal rate estimates and the consequential inference of neural interactions will mislead the coordination of the neural network. Hence, it will be necessary to analyze how discrepancies in rate estimates change the properties of DONAs or PDONAs and how the suitability of DONAs or PDONAs can be judged when using diverse datasets with varying levels of difficulty.
Second, timescales describing SNA and TSNA will vary according to types of neuron or task [22,38]. Furthermore, correlations between neurons are identified in accordance with the timescale of their pairwise effects [12]. Consequently, it is vital that the temporal precision of the rate estimates and the transient effects are chosen appropriately to allow confident interpretation of neural interactions. In the present study, we chose the temporal precisions heuristically so that the analysis model showed results that were more distinctive. Therefore, we cannot exclude the possibility that models with different temporal precisions might provide new information about neural interactions. Additionally, for real-time processing, such as in brain-machine interfaces or closed-loop brain stimulations, temporal precisions must be automatically regulated depending on the given neural data. Thus, for the general application of our method, it will be essential to identify techniques and criteria by which to define temporal precisions according to the characteristics of the neural data.
Third, we only compared our methods to cross-correlation coefficients. A number of other measures exist for quantifying interactions that have various uses related to the characteristics of neural data or the purposes of analysis [10][11][12][13][14]. To validate the effectiveness of DONA and PDONA as measures for inferring neural interaction, comparisons to other measures, including the degree of synchronous firing or the temporal correlation of neural pairs, will be required to identify the similarities and differences. Indeed, we should be open to the possibility that combinations of single-trial rate estimates and measures of correlations might provide more information about network structure.
We anticipate that this study will be extensible for further research into neural population dynamics. Our developed method can provide novel insights into types of neurons or neural pairs based on neural interactions, just as neurons are classified into certain types based on their functional roles or responses in specific epochs or sessions [43]. Here we applied our method only to spike data with binary values. For its application to continuous datasets, such as those containing EEG or fMRI data, further research will be required. In particular, we must be