State-Space Analysis of Time-Varying Higher-Order Spike Correlation for Multiple Neural Spike Train Data

Precise spike coordination between the spiking activities of multiple neurons is suggested as an indication of coordinated network activity in active cell assemblies. Spike correlation analysis aims to identify such cooperative network activity by detecting excess spike synchrony in simultaneously recorded multiple neural spike sequences. Cooperative activity is expected to organize dynamically during behavior and cognition; therefore currently available analysis techniques must be extended to enable the estimation of multiple time-varying spike interactions between neurons simultaneously. In particular, new methods must take advantage of the simultaneous observations of multiple neurons by addressing their higher-order dependencies, which cannot be revealed by pairwise analyses alone. In this paper, we develop a method for estimating time-varying spike interactions by means of a state-space analysis. Discretized parallel spike sequences are modeled as multi-variate binary processes using a log-linear model that provides a well-defined measure of higher-order spike correlation in an information geometry framework. We construct a recursive Bayesian filter/smoother for the extraction of spike interaction parameters. This method can simultaneously estimate the dynamic pairwise spike interactions of multiple single neurons, thereby extending the Ising/spin-glass model analysis of multiple neural spike train data to a nonstationary analysis. Furthermore, the method can estimate dynamic higher-order spike interactions. To validate the inclusion of the higher-order terms in the model, we construct an approximation method to assess the goodness-of-fit to spike data. In addition, we formulate a test method for the presence of higher-order spike correlation even in nonstationary spike data, e.g., data from awake behaving animals. The utility of the proposed methods is tested using simulated spike data with known underlying correlation dynamics. Finally, we apply the methods to neural spike data simultaneously recorded from the motor cortex of an awake monkey and demonstrate that the higher-order spike correlation organizes dynamically in relation to a behavioral demand.


Introduction
Precise spike coordination within the spiking activities of multiple single neurons is discussed as an indication of coordinated network activity in the form of cell assemblies [1] comprising neuronal information processing. Possible theoretical mechanisms and conditions for generating and maintaining such precise spike coordination have been proposed on the basis of neuronal network models [2][3][4]. The effect of synchronous spiking activities on downstream neurons has been theoretically investigated and it was demonstrated that these are more effective in generating output spikes [5]. Assembly activity was hypothesized to organize dynamically as a result of sensory input and/or in relation to behavioral context [6][7][8][9][10]. Supportive experimental evidence was provided by findings of the presence of excess spike synchrony occurring dynamically in relation to stimuli [11][12][13][14], behavior [14][15][16][17][18][19], or internal states such as memory retention, expectation, and attention [8,[20][21][22][23].
Over the years, various statistical tools have been developed to analyze the dependency between neurons, with continuous improvement in their applicability to neuronal experimental data (see [24][25][26] for recent reviews). The cross-correlogram [27] was the first analysis method for detecting the correlation between pairs of neurons and focused on the detection of stationary correlation. The joint-peri stimulus time histogram (JPSTH) introduced by [11,28] is an extension of the cross-correlogram that allows a time resolved analysis of the correlation dynamics between a pair of neurons. This method relates the joint spiking activity of two neurons to a trigger event, as was done in the peristimulus time histogram (PSTH) [29][30][31] for estimating the time dependent firing rate of a single neuron. The Unitary Event analysis method [25,32,33] further extended the correlation analysis to enable it to test the statistical dependencies between multiple, nonstationary spike sequences against a null hypothesis of full independence among neurons. Staude et al. developed a test method (CuBIC) that enables the detection of higher-order spike correlation by computing the cumulants of the bin-wise population spike counts [34,35].
In the last decade, other model-based methods have been developed that make it possible to capture the dependency among spike sequences by direct statistical modeling of the parallel spike sequences. Two related approaches based on a generalized linear framework are being extensively investigated. One models the spiking activities of single neurons as a continuous-time point process or as a discrete-time Bernoulli process. The point process intensities (instantaneous spike rates) or Bernoulli success probabilities of individual neurons are modeled in a generalized linear manner using a log link function or a logit link function, respectively [36][37][38]. The dependency among neurons is modeled by introducing coupling terms that incorporate the spike history of other observed neurons into the instantaneous spike rate [38][39][40]. Recent development in causality analysis for point process data [41] makes it possible to perform formal statistical significance tests of the causal interactions in these models. Typically, the models additionally include the covariate stimulus signals in order to investigate receptive field properties of neurons, i.e., the relations between neural spiking activities and the known covariate signals. However, they are not suitable for capturing instantaneous, synchronous spiking activities, which are likely to be induced by an unobserved external stimulus or a common input from an unobserved set of neurons. Recently, a model was proposed to dissociate instantaneous synchrony from the spike-history dependencies; it additionally includes a common, non-spike driven latent signal [42,43]. These models provide a concise description of the multiple neural spike train data by assuming independent spiking activities across neurons conditional on these explanatory variables. As a result, however, they do not aim to directly model the joint distribution of instantaneous spiking activities of multiple neurons.
In contrast, an alternative approach, which we will follow and extend in this paper, directly models the instantaneous, joint spiking activities by treating the neuronal system as an ensemble binary pattern generator. In this approach, parallel spike sequences are represented as binary events occurring in discretized time bins, and are modeled as a multivariate Bernoulli distribution using a multinomial logit link function. The dependencies among the binary units are modeled in the generalized linear framework by introducing undirected pairwise and higher-order interaction terms for instantaneous, synchronous spike events. This statistical model is referred to as the 'log-linear model' [44,45], or Ising/ spin-glass model if the model contains only lower-order interactions. The latter is also referred to as the maximum entropy model when these parameters are estimated under the maximum likelihood principle. In contrast to the former, biologically-inspired network-style modeling, the latter approach using a log-linear model was motivated by the computational theory of artificial neural networks originating from the stationary distribution of a Boltzmann machine [46,47], which is in turn the stochastic analogue of the Hopfield network model [48,49] for associative memory.
The merit of the log-linear model is its ability to provide a welldefined measure of spike correlation. While the cross-correlogram and JPSTH provide a measure of the marginal correlation of two neurons, these methods cannot distinguish direct pairwise correlations from correlations that are indirectly induced through other neurons. In contrast, a simultaneous pairwise analysis based on the log-linear model (an analogue of the Ising/spin-glass analysis in statistical mechanics) can sort out all of the pairdependencies of the observed neurons. A further merit of the loglinear model is that it can provide a measure of the 'pure' higherorder spike correlation, i.e., a state that can not be explained by lower-order interactions. Using the viewpoint of an information geometry framework [50], Amari et al. [44,45,51] demonstrated that the higher-order spike correlations can be extracted from the higher-order parameters of the log-linear model (a.k.a. the natural or canonical parameters). The strengths of these parameters are interpreted in relation to the lower-order parameters of the dual orthogonal coordinates (a.k.a. the expectation parameters). The information contained in the higher-order spike interactions of a particular log-linear model can be extracted by measuring the distance (e.g., the Kullback-Leibler divergence) between the higher-order model and its projection to a lower-order model space, i.e., a manifold spanned by the natural parameters whose higher-order interaction terms are fixed at zero [44,[52][53][54].
Recently, a log-linear model that considered only up to pairwise interactions (i.e., an Ising/spin-glass model) was proposed as a model for parallel spike sequences. Its adequateness was shown by the fact that the firing rates and pairwise interactions explained more than *90% of the data [55][56][57]. However, Roudi et al. [54] demonstrated that the small contribution of higher-order correlations found from their measure based on the Kullback-Leibler divergence could be an artifact caused by the small number of neurons analyzed. Other studies have reported that higher-order correlations are required to account for the dependencies between parallel spike sequences [58,59], or for stimulus encoding [53,60]. In [60], they reported the existence of triple-wise spike correlations in the spiking activity of the neurons in the visual cortex and showed their stimulus dependent changes. It should be noted, though, that these analyses assumed stationarity, both of the firing rates of individual neurons and of their spike correlations. This was possible because the authors restricted themselves to data recorded either from in vitro slices or from anesthetized animals. However, in order to assess the behavioral relevance of pairwise and higher-order spike correlations in awake behaving animals, it is necessary to appropriately correct for time-varying firing rates within an experimental trial and provide an algorithm that reliably estimates the time-varying spike correlations within multiple neurons.

Author Summary
Nearly half a century ago, the Canadian psychologist D. O. Hebb postulated the formation of assemblies of tightly connected cells in cortical recurrent networks because of changes in synaptic weight (Hebb's learning rule) by repetitive sensory stimulation of the network. Consequently, the activation of such an assembly for processing sensory or behavioral information is likely to be expressed by precisely coordinated spiking activities of the participating neurons. However, the available analysis techniques for multiple parallel neural spike data do not allow us to reveal the detailed structure of transiently active assemblies as indicated by their dynamical pairwise and higherorder spike correlations. Here, we construct a state-space model of dynamic spike interactions, and present a recursive Bayesian method that makes it possible to trace multiple neurons exhibiting such precisely coordinated spiking activities in a time-varying manner. We also formulate a hypothesis test of the underlying dynamic spike correlation, which enables us to detect the assemblies activated in association with behavioral events. Therefore, the proposed method can serve as a useful tool to test Hebb's cell assembly hypothesis.
We consider the presence of excess spike synchrony, in particular the excess synchrony explained by higher-order correlation, as an indicator of an active cell assembly. If some of the observed neurons are a subset of the neurons that comprise an assembly, they are likely to exhibit nearly completely synchronous spikes every time the assembly is activated. It may be that such spike patterns are not explained by mere pairwise correlations, but require higher-order correlations for explanation of their occurrence. One of the potential physiological mechanisms for higher-order correlated activity is a common input from a set of unobserved neurons to the assembly that includes the neurons under observation [25,[61][62][63]. Such higher-order activity is transient in nature and expresses a momentary snapshot of the neuronal dynamics. Thus, methods that are capable of evaluating time-varying, higher-order spike correlations are crucial to test the hypothesis that biological neuronal networks organize in dynamic cell assemblies for information processing. However, many of the current approaches based on the log-linear model [44,45,53,55,56,61,62,64,65] are not designed to capture their dynamics. Very recently two approaches were proposed for testing the presence of non-zero pairwise [66] and higher-order [67] correlations using a time-dependent formulation of a log-linear model. In contrast to these methods, the present paper aims to directly provide optimized estimates of the individual time-varying interactions with confidence intervals. This enables to identify short lasting, time-varying higher-order correlation and thus to relate them to behaviorally relevant time periods.
In this paper, we propose an approach to estimate the dynamic assembly activities from multiple neural spike train data using a 'state-space log-linear' framework. A state-space model offers a general framework for modeling time-dependent systems by representing its parameters (states) as a Markov process. Brown et al. [37] developed a recursive filtering algorithm for a point process observation model that is applicable to neural spike train data. Further, Smith and Brown [68] developed a paradigm for joint state-space and parameter estimation for point process observations using an expectation-maximization (EM) algorithm. Since then, the algorithm has been continuously improved and was successfully applied to experimental neuronal spike data from various systems [38,[69][70][71] (see [72] for a review). Here, we extend this framework, and construct a multivariate state-space model of multiple neural spike sequences by using the log-linear model to follow the dynamics of the higher-order spike interactions. Note that we assume for this analysis typical electrophysiological experiments in which multiple neural spike train data are repeatedly collected under identical experimental conditions ('trials'). Thus, with the proposed method, we deal with the within-trial nonstationarity of the spike data that is expected in the recordings from awake behaving animals. We assume, however, that dynamics of the spiking statistics within trials, such as time-varying spike rates and higher-order interactions, are identical across the multiple experimental trials (across-trial stationarity).
To validate the necessity of including higher-order interactions in the model, we provide a method for evaluating the goodness-offit of the state-space model to the observed parallel spike sequences using the Akaike information criterion [73]. We then formulate a hypothesis test for the presence of the latent, time-varying spike interaction parameters by combining the Bayesian model comparison method [74][75][76] with a surrogate method. The latter test method provides us with a tool to detect assemblies that are momentarily activated, e.g., in association with behavioral events. We test the utility of these methods by applying them to simulated parallel spike sequences with known dependencies. Finally, we apply the methods to spike data of three neurons simultaneously recorded from motor cortex of an awake monkey and demonstrate that a triple-wise spike correlation dynamically organizes in relation to a behavioral demand.
The preliminary results were presented in the proceedings of the IEEE ICASSP meeting in 2009 [77], as well as in conference abstracts (Shimazaki et al., Neuro08, SAND4, NIPS08WS, Cosyne09, and CNS09).

General formulation
Log-linear model of multiple neural spike sequences.
We consider an ensemble spike pattern of N neurons. The state of each neuron is represented by a binary random variable, X i (i~1, . . . ,N) where '1' denotes a spike occurrence and '0' denotes silence. An ensemble spike pattern of N neurons is represented by a vector, X~X 1 ,X 2 , . . . ,X N ð Þ , with the total number of possible spike patterns equal to 2 N . Let p(x), where x~x 1 ,x 2 , . . . ,x N ð Þ and x i~1 or 0 (i~1, . . . ,N), represent the joint probability mass function of the N-tuple binary random variables, X . Because of the constraint P p(x)~1, the probabilities of all the spike patterns are specified using 2 N {1 parameters. In information geometry [44,50], these 2 N {1 parameters are viewed as 'coordinates' that span a manifold composed of the set of all probability distributions, fp(x)g. In the following, we consider two coordinate systems.
The logarithm of the probability mass function can be expanded as where h~(h 1 ,h 2 , . . . ,h 12 ,h 13 , . . . ,h 1ÁÁÁN )'. In this study, a prime indicates the transposition operation to a vector or matrix. y(h)~{log p(f0, . . . ,0g) is a log normalization parameter to satisfy P p(x)~1. The log normalization parameter, y(h), is a function of h. Eq. 1 is referred to as the log-linear model. The parameters h i , h ij ,…, h 1ÁÁÁN are the natural parameters of an exponential family distribution and form the 'h-coordinates'. The h-coordinates play a central role in this paper.
The other coordinates, called g-coordinates, can be constructed using the following expectation parameters: These parameters represent the expected rates of joint spike occurrences among the neurons indexed by the subscripts. We denote the set of expectation parameters as a vector g~(g 1 ,g 2 , . . . ,g 12 ,g 13 , . . . ,g 1ÁÁÁN )'. Eqs. 1 and 2 can be compactly written using the following 'multi-indices' notation. Let V k be collections of all the k-element subsets of a set with N elements (i.e., a k-subset of N elements): We use the 'multi-indices' I[fV 1 , . . . ,V N g to represent the natural parameters and expectation parameters as h I and g I , respectively. Similarly, let us denote the interaction terms in Eq. 1 as Using Eqs. 1 and 2 can now be compactly written as log p(x)~P I h I f I x ð Þ{y(h) and g I~E ½f I (X), respectively.
The higher-order natural parameters in the log-linear model represent the strength of the higher-order spike interactions. Amari et al. [44,45,50,51] proved that the hand g-coordinates are dually 'orthogonal' coordinates and demonstrated that the natural parameters that are greater than or equal to the rth-order, , represent an excess or paucity of higher-order synchronous spikes in the fg To understand the potential influence of the higherorder natural parameters on the higher-order joint spike event rates, let us consider a log-linear model in which the parameters higher than or equal to the rth-order vanish: In the dual representation, the higher-order joint spike event rates, fg I g (I[ V r , . . . ,V N f g ), are chance coincidences expected from the lower-order joint spike event rates, fg I g (i[ V 1 , . . . ,V r{1 f g ). From this, it follows that the nonzero higher-order natural parameters that are greater than or equal to the rth-order of a full log-linear model represent the excess or paucity of the higher-order joint spike event rates, , in comparison to their chance rates expected from the lower-order joint spike event rates, fg I g (I[ V 1 , . . . ,V r{1 f g ). However, in this framework, the excess or scarce synchronous spike events of r subset neurons reflected in g I (I[V r ) are caused not only by the non-zero rth-order interactions fh I g (I[V r ), but also by all non-zero higher-order interactions fh I g (I[ V r , . . . ,V N f g ). Therefore, one cannot extract the influence of pure rth-order spike interactions (rvN) on the higher-order synchrony rates unless the parameters higher than the rth-orders vanish, h I~0 (I[ V rz1 , . . . ,V N f g ). To formally extract the rth-order spike interactions, let E r denote a log-linear model whose parameters higher than the rth-order are fixed at zero. By successively adding higher-order terms into the model, we can consider a set of log-linear models forming a hierarchical structure as E 1 5E 2 5 Á Á Á 5E N . Here, E r 5E rz1 denotes that E r is a submanifold of E rz1 in the h-coordinates. In E r , the set of non-zero rth-order natural parameters explains the excess or paucity of rthorder synchronous spike events, g I (I[V r ), in comparison to their chance occurrence rates expected from the lower-order marginal rates, . In a full model E N , the last parameter in the log-linear model, h 1ÁÁÁN , represents the pure Nth-order spike correlation among the observed N neurons. The information geometry theory developed by Amari and others [44,45,50,51] provides a framework for illustrating the duality and orthogonal relations between the hand g-coordinates and singling out the higher-order correlations from these coordinates using hierarchical models. In the Methods subsection 'Mathematical properties of log-linear model' we describe the known properties of the log-linear model utilized in this study, and give references [44,45,50,51] for further details on the information geometry approach to spike correlations.

0
, where each element In neurophysiological experiments, neuronal responses are repeatedly obtained under identical experimental conditions ('trials'). Thus, we here consider the parallel spike sequences repeatedly recorded simultaneously from N neurons over n trials and align them to a trigger event. We model these simultaneous spike sequences as time-varying, multivariate binary processes by assuming that these parallel spike sequences are discretized in time into T bins of bin-width D. Let X t,l~( X t,l 1 ,X t,l 2 , . . . ,X t,l N ) be the observed N-tuple binary variables in the t-th bin of the l-th trial, whereby a bin containing '1' values indicates that one or more spikes occurred in the bin whereas '0' values indicate that no spike occurred in the bin. We regard the observed spike pattern X t,l at time t as a sample (from a total of n trials) from a joint probability mass function. An efficient estimator of g t I is the observed spike synchrony rate defined for the t-th bin as for I[ V 1 , . . . ,V r f g . The synchrony rates up to the r-th order, y t~½ y t 1 , . . . ,y t 12 , . . . ,y t 1ÁÁÁr , . . .

0
, constitute a sufficient statistic for the log-linear model up to the r-th order interaction.
Using Eqs. 4 and 6, the likelihood of the observed parallel spike sequences is given as with y 1:T~y1 ,y 2 , . . . ,y T f gand h 1:T~h1 ,h 2 , . . . ,h T f g . Here, we assumed that the observed spike patterns are conditionally independent across bins given the time-dependent natural parameters, and that the samples across trials are independent.
Our prior assumption about the time-dependent natural parameters is expressed by the following state equation for t~2, . . . ,T. Matrix F (d|d matrix) contains the first order autoregressive parameters. j t (d|1 matrix) is a random vector drawn from a zero-mean multivariate normal distribution with covariance matrix Q (d|d matrix). The initial values obey h 1 *N m,S ð Þ, with m (d|1 matrix) being the mean and S (d|d matrix) the covariance matrix. The parameters F, Q, m, and S are called hyper-parameters. In the following, we denote the set of hyper-parameters by w: w~F,Q,m,S ½ . Given the likelihood (Eq. 7) and prior distribution (Eq. 8), we aim at obtaining the posterior density p h 1:T jy 1:T ,w ð Þ p y 1:T jh 1: using Bayes' theorem. The posterior density provides us with the most likely paths of the log-liner parameters given the spike data (maximum a posteriori, MAP, estimates) as well as the uncertainty in its estimation. The posterior density depends on the choice of the hyper-parameters, w. Here, the hyper-parameters, w, except for the covariance matrix S for the initial parameters, are optimized using the principle of maximizing the logarithm of the marginal likelihood (the denominator of Eq. 9, also referred to as the evidence), For non-Gaussian observation models such as Eq. 7, the exact calculation of Eq. 10 is difficult (but see the approximate formula in the Methods section). Instead, we use the EM algorithm [68,78] to efficiently combine the construction of the posterior density and the optimization of the hyper-parameters under the maximum likelihood principle. Using this algorithm, we iteratively construct a posterior density with the given hyper-parameters (E-step), and then use it to optimize the hyper-parameters (M-step). To obtain the posterior density (Eq. 9) in the E-step, we develop a nonlinear recursive Bayesian filter/smoother. The filter distribution is sequentially constructed by combining the prediction distribution for time t based on the state equation (Eq. 8) and the likelihood function for the spike data at time t, Eq. 4. Figure 1 illustrates the recursive filtering process in model subspace E r . In combination with a fixed-interval smoothing algorithm, we derive the smooth posterior density (Eq. 9). The time-dependent log-linear parameters (natural parameters) are estimated as MAP estimates of the optimized smooth posterior density. In the Methods section, please refer to the subsection on 'Bayesian estimation of dynamic spike interactions' for the derivation of the optimization method for hyper-parameters, along with the filtering/smoothing methods. We summarize a method for estimating the dynamic spike interactions in Table 1.
Application of state-space log-linear model to simulated spike data Estimation of time-varying pairwise spike interaction.
To demonstrate the utility of the developed methods for the analysis of dynamic spike correlations, we first consider a nonstationary pairwise spike correlation analysis. For this goal, we apply the state-space method to two examples of simulated spike data, with N~2 neurons. The dynamic spike correlation between two neurons can be analyzed by conceptually simpler histogram-based methods, e.g., a joint peri-stimulus time histogram (JPSTH). However, even for the pair-analysis, the proposed method can be advantageous in the following two aspects. First, the proposed method provides a credible interval (a Bayesian analogue of a confidence interval). Using the recursive Bayesian filtering/smoothing algorithm developed in the Methods section, we obtain the joint posterior density of the log-linear parameters (Eq. 9). The posterior density provides, not only the most likely path of the log-liner parameters (MAP estimates), but also the uncertainty in its estimation. The credible interval allows us to examine whether the pairwise spike correlation is statistically significant (but see the later section on ''Testing spike correlation in nonstationary spike data'' for the formal use of the joint posterior density for testing the existence of the spike correlation in behaviorally relevant time periods). Second, an EM algorithm developed in the proposed method optimizes the smoothness of the estimated dynamics of the pairwise correlation (i.e., optimization of the hyper-parameters, w, in the state equation, Eq. 8). By the automatic selection of the smoothness parameter, we can avoid the problem of spurious modulation in the estimated dynamic spike correlation caused by local noise, or excessive smoothing of the underlying modulation. Figure 2A displays an application of our state-space method to 2 parallel spike sequences, X t,l i (i~1,2, t~1, . . . ,T, and l~1, . . . ,n), which are correlated in a time-varying fashion. The data are generated as realizations from a time-dependent formulation of a full log-linear model of 2 neurons (Figure 2A left, repeated trials: n~50; duration: T~400 bins of width D). Here, the underlying model parameters, h t 1 , h t 2 , and h t 12 (Figure 2A right, dashed lines), are designed so that the individual spike rates are constant (g t 1~0 :0384 ½spike=D and g t 2~0 :0194 ½spike=D), while the spike correlation between the two neurons, h t 12 , varies in time, i.e., across bins (synchronous spike events caused by the timedependent correlation are marked as blue circles in Figure 2A left). While the bin-width, D, is an arbitrary value in this simulation analysis, the bin-width typically selected in spike correlation analyses is on the milli-second order. If the bin-width is D~1 ms, the individual spike rates of simulated neurons 1 and 2 are 38.4 Hz and 19.4 Hz, respectively. The correlation coefficient calculated from these parallel spike sequences is 0.0763. These values are within the realistic range of values obtained from experimentally recorded neuronal spike sequences. By applying the state-space method to the parallel spike sequences, we obtain the smooth posterior density of the log-linear parameters. The right panels in Figure 2A display the MAP estimates of the log-linear parameters (solid lines). The analysis reveals the time-varying pairwise interaction between the two neurons (Figure 2A right, bottom). The gray bands indicate 99% credible intervals from the posterior density, Eq. 9. We used marginal posterior densities, p(h t I jy 1:T ,w) (I[fV 1 ,V 2 g), to display the credible intervals for the individual loglinear parameters. The variances of the individual marginal densities were obtained from the diagonal of a covariance matrix of the smooth joint posterior density (Eq. 35). Figure 2B shows an application of the method to parallel spike sequences of time-varying spike rates. Here, the underlying model parameters are constructed so that the two parallel spike sequences are independent (h t 12~0 for t~1, . . . ,T), while the individual spike rates vary in time (i.e., across the bins). The observation of synchronous spike events ( Figure 2B left, blue circles) confirms that chance spike coincidences frequently and trivially occur at higher spike rates. The analysis based on our state-space method reveals that virtually no spike correlation exists between the two neurons, despite the presence of time-varying rates of synchronous spike events ( Figure 2B right, bottom).

Simultaneous estimation of time-varying pairwise spike
interactions. In this subsection, we extend the pairwise correlation analysis of 2 neurons to the simultaneous analysis of multiple pairwise interactions in the parallel spike sequences obtained from more than 2 neurons. Figure 3 demonstrates an application of our method to simulated spike sequences of N~8 neurons. As an underlying model, we construct a time-dependent log-linear model of 8 neurons with time-varying rates and pairwise interactions (h t I , I[V 1 ,V 2 , duration: T~500 bins). The higher-order log-linear parameters are set to zero, i.e., no higher-order interactions are included in the model. Figure 3A displays snapshots of the dynamics of the parameters of the individual spike rates, g t I (I[V 1 ), and pairwise interactions, h t I (I[V 2 ), at t~100,200,300, 400 and 500. Figure 3B shows the parallel spike sequences (50 out of 200 trials are displayed) simulated on the basis of this model. The spikes involved in the pairwise, synchronous spike events between any two of the neurons (in total: 28 pairs) are superimposed and marked with blue circles. Figure 3C displays snapshots of the simultaneous MAP estimates of the pairwise interactions, h t I (I[V 2 ), of a log-linear model of 8 neurons applied to the parallel spike train data. In addition, the spike rates were estimated from the dual coordinates, i.e., g t I (I[V 1 ). The results demonstrate that the simultaneous estimation of time-varying, multiple pairwise interactions can be carried out by using a statespace log-linear model with up to pairwise interaction terms.
Estimation of time-varying triple-wise spike interaction.
Another important aspect of the proposed method is its ability to estimate time-varying higher-order spike interactions that cannot be revealed by a pairwise analysis. To demonstrate this, we apply the state-space log-linear model to N~3 parallel spike sequences by considering up to a triple-wise interaction (i.e., the full log-linear model). Spike data ( Figure 4A) are generated by a time-dependent log-linear model ( Figure 4C, dashed lines) repeatedly in n~100 Table 1. Method for estimating dynamic spike interactions.

I Preprocessing parallel spike data
(1) Align parallel spike sequences from N neurons at the onset of external clock that repeated n times.
(3) Select r~1,2, . . . ,N ð Þ , the order of interactions included in the model. At each bin, compute the joint spike event rates up to the rth order y t , using Eq. 6. (3) M-step: Optimize the hyper-parameters.

II
Update the hyper-parameters, F and Q, using Eqs. 38 and 39, and m~h 1jT .
(4) Repeat (2) and (3)  trials. Figure 4C displays the MAP estimates (solid lines) of the loglinear parameters from the data shown in Figure 4A. Here, nonzero parameter h t 123 represents a triple-wise spike correlation, i.e., excess synchronous spikes across the three neurons or absence of such synchrony compared to the expectation if assuming pairwise correlations. The gray band is the 99% credible interval from the marginal posterior density, p h t I jy 1: The credible interval of the higher-order (triple-wise) log-linear parameter in the bottom panel of Figure 4C appears to be larger than those of the lower-order parameters. In general, the observed frequency of simultaneous spike occurrences decreases as the number of neurons that join the synchronous spiking activities increases (note that the marginal joint spike occurrence rate, Eq. 2, is a non-increasing function with respect to the order of interaction, i.e., g I §g J if the elements of I are included in J). Thus, the estimation variance typically increases for the higherorder parameters, as seen in the bottom panel of Figure 4C. Related to the above, because of the paucity of samples for higher-order joint spike events, the automatic smoothness optimization method selects hyper-parameters that make the trajectories of the higher-order log-linear parameters stiff in order to avoid statistical fluctuation caused by a local noise structure. Given the limited number of trials available for data analyses, these observations show the necessity of a method to validate inclusion of the higherorder interaction terms in the model.
Additionally, we observe that in the later period of spike data (300-500 bins), the dynamics of the estimated triple-wise spike interaction do not follow the underlying trajectory faithfully: The underlying trajectory falls on outside the 99% credible interval. Similar results are sometimes observed when an autoregressive parameter, F, in a state model is optimized (Eq. 8). In contrast, when we replace the autoregressive parameter with an identity matrix (i.e., F~I, where I is the identity matrix), the credible intervals become larger. Therefore, such observations do not typically occur. Thus, we also need a method for validating the inclusion of the autoregressive parameter in the state model using an objective criterion. Detailed analyses of these topics will be given in the next section using the example of 3 simulated neurons displayed in Figure 4. In the above example, the state-space loglinear model with an optimized F provides a better overall fits to the spike data than the model using F~I, despite an inaccurate representation in part of its estimation. However, for the purpose of testing the spike correlation in a particular period of spike data, we recommend using an identity matrix as an autoregressive parameter, i.e. F~I, in the state model.

Selection of state-space log-linear model
For a given number of neurons, N, we can construct state-space log-linear models that contain up to the rth-order interactions (r~1, . . . ,N). While the inclusion of increasingly higher-order interaction terms in the model improves its accuracy when describing the probabilities of 2 N spike patterns, the estimation of the higher-order log-linear parameters of the model may suffer from large statistical fluctuations caused by the paucity of synchronous spikes in the data, leading to an erroneous estimation of such parameters. This problem is known as 'over-fitting' the model to the data. An over-fitted model explains the observed data, but loses its predictive ability for unseen data (e.g., spike sequences in a new trial under the same experimental conditions). In this case, the exclusion of higher-order parameters from the model may better explain the unseen data even if an underlying spike generation process contains higher-order interactions. The model that has this predictive ability by optimally resolving the balance between goodness-of-fit to the observed data and the model simplicity is obtained by maximizing the cross-validated likelihood or minimizing the so-called information criterion. In this section, we select a state-space model that minimizes the Akaike information criterion (AIC) [73], which is given as The first term is the log marginal likelihood, as in Eq. 10. The second term that includes k is a penalization term. The AIC uses the number of free parameters in the marginal model (i.e., the number of free parameters in w) for k. Please see in the Methods subsection 'Selection of state-space model by information criteria' for an approximation method to compute the marginal likelihood. Selecting a model that minimizes the AIC is expected to be equivalent to selecting a model that minimizes the expected (or average) distance between the estimated model and unknown underlying distribution that generated the data, where the 'distance' measure used is the Kullback-Leibler (KL) divergence. The expectation of the KL divergence is called the KL risk function. Selection from hierarchical models. Here, we examine the validity of including higher-order interaction terms in the model by using the AIC. We apply the model selection method to the spike train data of N~3 simulated neurons. The data are generated by the time-varying, full log-linear model that contains a non-zero triple-wise interaction terms shown in Figure 4C (dashed lines). The AICs are computed for hierarchical state-space loglinear models, i.e., for models of interaction orders up to r~1,2,3.
To test the influence of the data sample size on the model selection, we vary the number of trials, n, used to fit the hierarchical log-linear models. The results are shown in Table 2 for n~2, 5,20,100,200. For a small number of trials (n~5), a model without any interaction structure (r~1) is selected. For larger numbers of trials, models with larger interaction orders are selected. For n~100 and 200, the full log-linear model (r~3) is selected.
Below, we examine whether the AIC selected a model that minimizes the KL risk function by directly computing its approximation using the known underlying model parameters. First, Table 3 shows how often a specific order is selected by the AIC by repeatedly applying the method to different samples generated from the same underlying log-linear parameters ( Figure 4C, dashed lines). We examine two examples: One in which a sample is composed of n~5 trials (left) and the other of n~100 trials (right). We repeatedly compute the AICs of statespace models of different orders (r~1,2,3) applied to 100 data realizations (of the respective number of trials). We then count how often a model of order r is selected by minimizing the AIC. For comparison, the table includes the outcomes from other criteria such as the Bayesian information criterion (BIC) [79,80] and the predictive divergence for indirect observation models (PDIO) [81], which are suggested for models containing latent variables. Please see the Methods section for the details of these criteria. Next, Table 4 displays the most frequently selected model (from r = 1,2,3) by the various information criteria when they are applied to 100 data realizations as a function of the number of trials in each data set, n~2,5,20,100,200 (see Table 3 for the outcomes of n~5 and n~100).
We compare the outcomes of the AIC and two other information criteria with the model order that minimizes the KL risk function (KL-risk). For this goal, we include in Table 4 the model order that minimizes the KL-risk between the underlying log-linear model ( Figure 4C, dashed lines) and estimates of the r th-order model. In addition to the KL-risk, we calculate the mean squared error (MSE) between the underlying model parameters and corresponding estimates. Please see in the caption of Table 4 how we compute the KL-risk and MSE for this analysis. We find that the KL-risk and MSE select the same model, except for the case of n~5 trials. In comparison to the KL-risk, the BIC tends to select models with an excessively higher-order of interaction (overfitting) for a small number of trials (n~2) and tends to choose lower-order models for n~20 and 100. The PDIO mostly selects lower-order models. In contrast, the AIC follows the selection of the KL-risk minimization principle, except for n~20, where it shows a conservative choice compared to the KL-risk.
We repeat the same analysis for spike data generated from an underlying model that contains up to pairwise interactions, but does not contain the triple-wise interaction term. The purpose of this analysis is to show that the methods do not select models with excessively higher orders of interaction than those actually contained in the data. To construct such an underlying model, we project the full model (shown in Figure 4C) onto the subspace of a pairwise loglinear model, E 2 . The projection model does not contain any triplewise correlation, while the 1st and 2nd order expectation parameters (g I for I[fV 1 ,V 2 g) are the same as those of the full model that was used to generate the data in the analysis of Tables 2-4. Table 5 displays the most frequently selected model orders by the AIC, BIC, PDIO, along with the selections by the KL-risk and MSE. We find that the pairwise model is the most frequently selected model under all of the criteria for the samples with a large number of trials (n~100,200). Under this condition, only the AIC among the other data-driven methods follows the KL-risk selection. These results lead us to the conclusion that the AIC is a reliable measure to assess the goodness of fit of the state-space log-linear model.
Selection of state transition model. In addition to validating the inclusion of the order of spike interactions, we examine state models (Eq. 8) with different conditions for the hyper-parameters by the AIC, using as an example the same spike train data of 3 neurons with n~100 trials (displayed in Figure 4A). The tested state models are (i) a time-independent model in which ñ5 ñ20 n~100 n~200 The state-space log-linear models of different orders (r~1,2,3) are applied to samples of the N~3 spike sequences of n repeated trials generated from a time-dependent full log-linear model (indicated by the dashed lines in Figure 4C). Three data-driven information criteria, AIC, BIC, and PDIO, are computed for the fitted state-space models of the different orders, r~1,2,3.
The count for the model order r that minimizes the respective criteria is determined by repeating the process for 100 repetitions as in Table 3. The most frequently selected model order, r, is displayed for each of the information criteria and for the different numbers of trials, n. For comparison, we also show the order of interactions that minimizes the KL risk function (KL-risk) and mean squared error (MSE). We approximated the KL-risk and MSE as follows. At each bin, we compute the KL-divergence (Eq. 21), between a full underlying loglinear model of N~3 neurons and the estimated log-linear model whose parameters are given by the MAP estimates of the r th-order model. The total sum of the all KL-divergences from T~500 bins is used as the distance between the two (time-dependent) models: i.e., KL~P T t~1 D½q(h t ),p(h tjT ), where the function D½q,p is given in Eq. 21. h t represents the underlying loglinear parameters used to generate the data. h tjT is its estimate from one sample composed of n trials. The parameters higher than the r th-order that are not included in the model are set to zero. The KL-risk function is estimated as the average of the KL-divergences of 100 realizations of the spike data, each composed of n trials. To obtain the MSE, we first computed the sum of the squared errors (SE) : SE~P T t~1 jjh t {h tjT jj 2 , using one sample composed of n trials. The MSE is then estimated as the average of the SEs over 100 samples. doi:10.1371/journal.pcbi.1002385.t004  The state-space log-linear models containing interactions up to the rth-order (r~1,2,3) are applied to the data from N~3 simultaneous spike sequences. The spike data is generated from a time-dependent full log-linear model (see dashed lines in Figure 4C) repeatedly for n trials; either n~5 (left) or n~100 (right) trials. For this data set, we compute three information criteria (AIC, BIC, and PDIO) and find the order of the model that minimizes these information criteria. We repeated the selection of the model order 100 times, using each of the criteria and using the spike data that contains respective number of trials (n~5 or n~100). The count of the order of spike interactions that minimizes the applied information criteria is increased by 1, and finally expressed as a frequency. The asterisk marks the most frequently selected model. doi:10.1371/journal.pcbi.1002385.t003 the hyper-parameters are fixed as Q~0 and F~I, where I is an identity matrix; (ii) a random walk model (F~I), in which only the covariance matrix, Q, is optimized via the EM algorithm; and (iii) a 1st-order autoregressive model, in which Q and F are both optimized. With these settings, the state-space log-linear model of case (i) becomes a stationary log-linear model, which has been frequently employed in spike data analyses, often in the form of a maximum entropy model [55][56][57]82]. In contrast, the state-space models of cases (ii) and (iii) are nonstationary because the joint distribution of the spike pattern changes in time as a result of the time-varying log-linear parameters. The interaction parameters themselves are assumed to be nonstationary in the state model in case (ii), and can be either stationary or nonstationary in the state model in case (iii). To see this, let l i be the eigenvalues of F, i.e., the solutions for det(F{lI)~0. For a stationary process, the eigenvalues have to satisfy l i j jv1, otherwise the process is nonstationary. This excludes the case F~I.
The AICs of the full model with the various state-equations become progressively smaller for increasingly complex state-space models (38525 for case (i), 36263 for case (ii), and 36231 for case (iii)). The fact that conditions (ii) and (iii) are better fitted to the data than condition (i) confirms that the proposed method of timevarying spike correlation analysis performs better than an analysis based on a stationary log-linear model for this data. In addition, the fact that condition (iii) better fits the data than condition (ii) supports the use of autoregressive models (note: The estimation in Figure 4B is done with state model (iii)). We find that the fitted F (to the data in Figure 4A) yields an eigenvalue larger than 1, indicating that the underlying state process (T~500 bins) is modeled as a nonstationary process. We consider that the nonstationary state model is selected because of the relatively short observation period, during which a nonstationary trend for the parameters can appear even if the state process is stationary in the long run.
Test for presence of spike correlation in nonstationary spike sequences One of the goals of a time resolved analysis of spike correlation is to discover dynamical changes in the correlated activities of neurons that reflect the behavior of an animal. This implies the necessity of dealing with the within-trial nonstationarity that is typically present in the data from awake behaving animals. However, we know from other correlation analysis approaches that, if not well corrected for, nonstationary spike data bear the potential danger of generating false outcomes [25,83]. Here we deal with the within-trial nonstationarity of the data by using the state-space log-linear model while assuming identical dynamic spiking statistics across trials (across-trial stationarity). In order to correctly detect the time-varying correlation structure within trials, we apply to the state-space log-linear model a Bayesian model comparison method based on the Bayes factor (BF) [74][75][76], and combine it with a surrogate approach. The BF is a likelihood ratio for two different hypothetical models of latent signals, e.g., in our application, different underlying spike correlation structures. Using the BF, we determine which of the two spike correlation models the spike data supports. By computing the BF for a particular task period in a behavioral experiment, we can test whether the assumed correlation structure appears in association with the timing of the animal's behavior. In the following, we denote a specific task period of interest by the time period ½a,b.
In this study, the BF, B 12 (y a:b ), is defined as the ratio of the marginal likelihoods of the observed spike patterns, y a:b , in the time period ½a,b under different models, M 1 or M 2 , assumed for the hidden state parameters, By successively conditioning the past, the BF is computed by the multiplication of the bin-by-bin one-step BF given at time t as Here, the bin-by-bin BF at time t, B 12 (y t jy a:t{1 ), can be calculated as (see the Methods subsection, 'Bayesian model comparison method for detecting spike correlation'), where S i is the space of the interaction parameters, h t , for the model, M i (i~1,2). In Eq. 13, p h t jy 1:t ,w ð Þis the filter density and p h t jy 1:t{1 ,w ð Þis called the one-step prediction density, both of which are obtained in the Bayesian recursive filtering algorithm developed in the Methods section (cf. Eqs. 25,26 and Eqs. 31,32). Therefore, the bin-by-bin BF at time t, B 12 (y t jy 1:t{1 ), is the ratio of the odds (of opposing models) found by observing the spike train data up to time t (filter odds, the numerator in Eq. 13) to the odds predicted from y 1:t{1 without observing the data at time t (prediction odds, the denominator in Eq. 13). Thus, an unexpected synchronous spike pattern that significantly updates the filter odds for the interaction parameters from their predicted odds gives rise to a large absolute value for the BF. Because the posterior densities are approximated as a multivariate normal distribution in our filtering algorithm, the BF at time t can be easily computed by using normal distribution functions. Please see the subsection, 'Bayesian model comparison method for detecting spike correlation', in the Methods section for the derivation of Eq. 13 and detailed analysis of the BF.
The BF becomes larger than 1 if the data, y a:b , support model M 1 as opposed to M 2 as an underlying spike correlation structure and becomes smaller than 1 if the data support model M 2 as opposed to M 1 . Alternatively, it is possible to use the logarithm of the BF, known as the 'weight of evidence' [75] which becomes positive if the data support model M 1 as opposed to model M 2 Similar to Table 4, the state-space log-linear models (r~1,2,3) are repeatedly applied to 100 samples of N~3 spike sequences with n trials. In contrast to Table 4, here a sample of spike sequences is generated from a pairwise (r~2) time-dependent log-liner model of N~3 neurons. This underlying model is a projection of the full log-liner model examined in Table 4 (dashed lines in Figure 4C) to the pairwise model space. The frequencies of the r th order model that minimizes AIC, BIC, and PDIO are counted by repeatedly applying models with different orders (r~1,2,3) to 100 samples (each with n trials). The most frequently selected models are displayed for the different criteria and for different numbers of trials, n, as in and negative in the opposite situation. Below, we display the results for the BF in bit units (logarithm of the BF to base 2), i.e., the weight of evidence. By sequentially computing the bin-by-bin BF, we can obtain the weight of evidence in a period ½a,b as the summation of the local weight of evidence: log 2 B 12 (y 1:T )P b t~a log 2 B 12 (y t jy 1:t{1 ). An intuitive interpretation of the BF values is provided in the literature [74,76]. For example, in [76], a BF (weight of evidence) from 1.6 to 4.3 bit was interpreted as 'positive' evidence in favor of M 1 against M 2 . Similarly, a BF from 4.3 to 7.2 bit was interpreted as 'strong' evidence, and a BF larger than 7.2 bit was found to be 'very strong' evidence in favor of M 1 against M 2 . While the classical guidelines are useful in practical situations, they are defined subjectively. Thus, in this study, in order to objectively analyze the observed value of the BF, we combine the Bayesian model comparison method with a surrogate approach. In this surrogate method, we test the significance of the observed BF for the tested spike interactions by comparing it with the surrogate BFs computed from the null-data generated by destroying only the target spike interactions while the other structures such as the time-varying spike-rates and lower-order spike interactions are kept intact.
The BF in a behaviorally relevant sub-interval ½a,b can be computed from the optimized state-space log-linear model fitted to the entire spike train data in ½0,T. Here, for the purpose of testing spike correlation in the sub-interval, we recommend to use F~I in the state model because the autoregressive parameters are optimized for entire spike train data, which are not necessarily optimal for the sub-interval. Similarly, a typical trial-based experiment is characterized by discrete behavioral or behaviorally relevant events, e.g. movement onset after a go signal or a cue signal for trial start, etc. Thus, on top of the expected smooth timevarying change in the spike-rate and spike-correlation, sudden transitions may be expected in their temporal trajectories. Because we use time-independent smoothing parameters (i.e., hyper parameter Q in Eq. 8) that were optimized to entire data (see the EM algorithm in the Methods section), such abrupt changes may not be captured very well. This may cause a false detection or failure in the detection of the spike correlation at the edge of a task period. For such data, we suggest applying the Bayesian model comparison method to state-space models which are independently fitted to each of the task periods (or relatively smooth subintervals within each task period).
Detecting triple-wise spike correlation in simulated nonstationary spike sequences. In this subsection, we examine a method to test for the presence of a higher-order (triple-wise) spike correlation using simulated spike train data with known spike interaction dynamics. In this context, we again use the simulated spike data of N~3 neurons (of length 500 bins with n~100 trials) shown in Figure 4A, which was generated by the time-varying model, as shown in Figure 4C (dashed lines). In the following, we assume for this data an underlying experimental protocol that can be segmented into behaviorally-relevant time periods, I-IV, as indicated in Figure 5A, as for example a period before the trial starts, a preparatory period, a movement period, etc. (e.g., see [8,19] for a typical experimental protocol).
In each of the time periods, we compare opposing models for the hidden log-linear parameters of N~3 neurons, using a full model (r~3). In one hypothetical model, M 1 , we assume that the triplewise interaction term is positive, S 1 : h t 123 w0. The time (bin) index t expresses that the BF is computed under the same models for all of the time steps, t, in the respective time period. In the other model M 2 , we assume that the triple-wise interaction term is smaller than or equal to zero, S 2 : h t 123 ƒ0. Neither model makes specific assumptions about the first and second order log-linear parameters; thus, they are allowed to be real numbers. These parameters are integrated out in Eq. 13. In this simulation study, we independently applied a state-space log-linear model to each of the four periods and computed the respective BFs. Figure 5B displays the ground truth of the time-dependent triple-wise interaction parameter, h t 123 , of the full log-linear model. The gray areas indicate periods where the triple-wise spike interaction, h t 123 , is positive: Model M 1 is true in this period. In the remaining periods (white), the model of a negative triple-wise spike interaction, M 2 , is true.
We compare the observed BF for the tested correlation models to the BFs computed from surrogate data sets resampled under a null hypothesis of no tested order of spike interaction. Specifically, to test the existence of a triple-wise spike correlation, the surrogate BFs are computed from resampled spikes generated from a model containing no triple-wise spike interaction. For this, we first apply a pairwise state-space log-linear model (r~2) to the observed spike data to derive estimates of the time-varying spike rates and dynamic pair-interactions. Then we generate 1000 surrogate samples of N~3 parallel spike sequences, where each sample is composed of n~100 trials, using the fitted pairwise state-space log-linear model. Thus, the resampled spike sequences exhibit the same time-varying spike rates and pairwise correlations as the original data, but the triplet coincidences occur on a chance level, as expected by the individual spike rates and pairwise correlations among the neurons. The surrogate BFs of a triple-wise spike correlation (using S 1 : h t 123 w0, S 2 : h t 123 ƒ0) are computed for each surrogate data set by applying the full model (r~3). Figure 5C shows the observed BF (red vertical line and red triangle) and the cumulative distribution function (CDF) of the surrogate BFs (solid black lines) for each of the four time periods, I-IV. In each of the graphs, the gray area indicates the 95% confidence interval derived from the distribution of 1000 respective surrogate BFs. If the observed BF falls outside the confidence interval, the null hypothesis that no tested order of spike interaction exists is rejected. Further, by the two-tailed test, if the observed BF falls outside the confidence interval into the red area of the distribution on the right (i.e., significantly positive BF values), the data supports model M 1 as an underlying correlation structure, whereas if it falls into the blue area on the left, we conclude that it is a significantly negative BF, i.e., M 2 is supported. If the observed BF falls into the 95% confidence interval, the null hypothesis that no tested order of spike correlation exists is not rejected. We first look at the results for periods II and IV in Figure 5C. In these periods, the null hypothesis of no triple-wise spike correlation is rejected, and the observed BF correctly suggests a positive triple-wise correlation model (M 1 ), reflecting the fact that the underlying triple-wise interaction term is positive throughout the two periods (see periods II and IV in Figure 5B). In contrast, in periods I and III, where a negative triple-wise correlation model (M 2 ) is true, the null hypothesis of no triple-wise correlation cannot be rejected.
It is possible that the synchronous spike events observed among multiple neurons are detected as simultaneous increases in pair interaction terms using the log-linear model of the second order (r~2), without having to take higher-order spike interactions into consideration. Paradoxically, this issue becomes relevant when a triple-wise spike correlation is a dominant factor in the generation of synchronous spike events because the projection of a full model (r~3) that contains a positive triple-wise interaction term to a pairwise model induces simultaneous increases in the pair interaction terms in that projected model (r~2). In such a case, an analysis of the higher-order spike correlations might just add redundant information about an animal's behavior to lower-order analyses. Thus, we repeat an analysis similar to the one above under this alternative hypothesis, using a model that contains up to pairwise interactions (r~2). In this test, the BF (Eq. 12) is computed using the following assumed models on the hidden loglinear parameters: M 1 , all three neurons simultaneously exhibit positive pairwise interactions (S 1 : h t I w0 for all I[V 2 ); and M 2 , at least one of the pairwise interactions is not positive (S 2 : h t I ƒ0 for at least one I[V 2 ). The first order log-linear parameters are again assumed to be real values. Surrogate data sets are obtained by resampling spike sequences from the first order state-space loglinear model (r~1) fitted to the data. The surrogate BFs for the model of simultaneously positive pairwise interactions, as opposed to its complementary model, are computed for each of the surrogate data sets by applying a pairwise model (r~2). Figure 5D displays the ground truth of the dynamics of the pairwise interaction terms, h t I (I[V 2 ), of the pairwise log-linear model. These are obtained by projecting the full underlying model ( Figure 4C dashed lines) to the probability space of the secondorder model (E 2 ). Note that in period II, the pair-interaction terms simultaneously increase because the effect of the triple-wise interaction in E 3 is projected to this sub-space model (E 2 ) (see Figure 4C, middle panel, for comparison, and observe that there is no such increase in the pair-interaction terms of the full log-linear model). The gray areas indicate the period where model M 1 (simultaneously positive pairwise interactions) is true. Figure 5E displays the results for testing the simultaneous increases in the pair interactions. We again look at periods II and IV. In period II, the null hypothesis of independent spike sequences is rejected, and model M 1 is supported, reflecting the fact that the projected pairwise interactions are positive for most of this time period. In contrast, the same null hypothesis cannot be rejected in period IV because neither model M 1 nor M 2 alone can support the data in this period: Models M 1 and M 2 are true in the first and second halves of period IV, respectively (see Figure 5D).  Figure 5F summarizes the results for periods I-IV using bars in corresponding colors. The results of our tests reflect the dynamics of the underlying parameters in Figures 5B and D, and we find that in period IV, only the test for a triple-wise spike correlation detected the presence of interactions among all three neurons in the data. This is because, in this time period, the dynamics of the underlying triple-wise interaction correlate well with the assumed behavioral time table, whereas the dynamics of the simultaneous pairwise interactions do not. In summary, it is possible that the higher-order analysis allows us to discover the correlated activities of multiple neurons associated with behavioral events, which may not be revealed by pairwise analyses.
Detecting triple-wise spike correlation in neural spike data from behaving monkey. Finally, we demonstrate the application of our method on simultaneous spike recordings from the primary motor cortex (MI) of an awake, behaving monkey. The data were recorded by Alexa Riehle and her colleagues [8] to test the hypothesis that neuronal cooperativity is involved in the planning of motor actions. Therefore, Riehle et al. designed a behavioral task in which different durations of preparation intervals were provided to the monkeys before they had to perform an arm movement and touch a target on a screen. The detailed time table of the task is as follows (cf. Figure 6A). A trial started by a signal indicating to the monkey that he may initiate the task by pressing a button. After initiating the trial, the monkey had to wait for 1000 ms until a preparatory signal (PS) was delivered, indicating to the monkey that now the preparation interval started. After appearance of a response signal (RS) at 600, 900, 1200, or 1500 ms (randomly selected with equal probability), the monkey had to move his arm and touch the target position. The reaction times of the monkey showed a dependence on the duration of the preparatory period (PP): The longer the PPs the shorter the reaction times [8,19,84]. The conditional probability for the occurrence of the RS increases during the longer PP, in particular at the times when the RS may occur. The hypothesis for the underlying network activity realizing reduced reaction times with longer PPs was that the system is better prepared for the requested movement by increasingly enhancing the synchrony in the relevant network to facilitate the response. Indeed, Riehle et al. showed that cortical neurons modulate the degree of excess synchronous spiking activities independently from their individual firing rates in conjunction with the occurrence of the expected RSs [8,19,84]. An open question that could not be conclusively answered in that study was if larger groups of neurons than pairs coordinate their activity in relation to the expected events. Here, we approach this question by applying our newly developed method and will indeed demonstrate the existence of higher-order (triple-wise) spike interaction in relation to the behavioral task.
To that end, we analyze data which were in part analyzed by the Unitary Events analysis method (neuron id 2 and 3; shown in Figure 2 in [8]). Figure 6B shows this particular data set that consists of three simultaneously recorded neurons (neuron id 2, 3, and 5; Here we denote them as neuron 1, 2, and 3) observed during trials (n~36) of the longest preparatory period of 1500 ms. As in Riehle et al. (Figure 2 in [8]) we align the trials at the PS, and analyze the data for the interaction parameters as a function of time from 200 ms before the PS and 1800 ms after the PS. This time segment is composed of three behavioral epochs: an interval before the PS (200 ms), the PP with three expected signals (ESs, at 600, 900 and 1200 ms), and a period after the RS at 1500 ms including the reaction time and partly the movement (reaction and movement time, RT-MT) (see Figure 6A). The average firing rates of the neurons during the PP are 29.4, 12.9, and 41.9 Hz for neurons 1, 2, and 3 (neuron id 2, 3, and 5) respectively. We constructed binary sequences, X t,l , from the spike times of the neurons by binning the data using a bin-width of D~3 ms. Figure 6C displays the occurrence rates for the individual spiking activities (upper panel) and the pair and triple synchronous spike events (middle and bottom panels) detected in the bins with a width of D~3 ms.
We apply a full log-linear model of N~3 neurons with up to a triple-wise interaction term to the binary spike data shown in Figure 6C. Figure 6D displays the estimated dynamics of the interaction parameters in the log-linear model by using the method summarized in Table 1. The thick bold lines indicate the MAP estimates, i.e., the most probable paths, of the interaction parameters. The gray bands and thin solid lines mark the 95% credible interval (the Bayesian analogue of the confidence interval) computed from a marginal posterior density. The last parameter of the log-linear model, h t 123 , indicates the triple-wise interaction among the three neurons in the MI. Strong positive (or negative) value of this term indicates that the three neurons are triple-wise correlated, i.e. dependent in a manner that cannot be explained by pair correlations of the neurons. Specifically, a positive triple-wise interaction value means that the occurrence of the synchronous events among the three neurons is more frequent than the chance level expected from the observed individual rates and pairwise spike correlations among them.
The analysis of the MI neurons using the state-space log-linear model reveals that the triple-wise interaction, h t 123 , during the PP gradually increases, with additional local peaks at the ESs that also increase in heights towards the end of the PP ( Figure 6D, bottom panel). This result is consistent with the results found by Riehle and collaborators for a subset of the neurons of the same data set ( Figure 2 in [8]). In this previous study neuron 1 and 2 (neuron id 2, 3) were analyzed for excess spike synchrony using the Unitary Events analysis [32,33]. The analysis revealed that the two neurons exhibit a modulation of significant excess spike synchrony, with peaks at the ESs (despite the first) and at the RS. [8,32,84]. That observation was interpreted as evidence that the neurons cooperate to prepare for motor action and facilitate the efficiency for movement execution [8,84]. Similarly, occurrences of synchronous spike events of the three neurons roughly coincide with the expected events (ESs) ( Figure 6C, bottom panel); accordingly, our result shows that the triple-wise interaction, h t 123 , is also locked to the ESs, however decreases at RS. Note that the triple-wise interaction is not only determined by the frequency of synchronous spike events of all three neurons but is also determined by frequencies of other observed spike patterns across trials. However, for typical neural spike train data with low spike rates and low rates of synchronous spiking of pairs of neurons, the occurrence of synchronous spike events of three neurons can significantly increase the triple-wise interaction. In the subsection 'Bayesian model comparison method for detecting spike correlation' in the Methods section we explore the contribution of different spike patterns in different spiking scenarios to the evidence of a triple-wise spike correlation in simulated data. Figure 6E displays modulation of the interaction term, h t 123 , using the binary data constructed with bins of smaller (D~2 ms, upper panel) and larger (D~5 ms, lower panel) widths. The top of each panel shows the timing of the synchronous spike events of all three neurons. The sample sizes for synchronous spike events across all three neurons for the 2 ms bin-width are greatly reduced from those observed for a bin-width of 3 ms. Thus estimated dynamics is much less structured: The data prevents us from detecting existence of a triple-wise spike correlation using the proposed statistical test (see below for the application of the test method). With larger bins of a width of 5 ms, the precise locking of the  [8,84]. During the preparatory period (PP, 1500 ms) that starts with the preparatory signal (PS), the presentation of the response signal (RS) was expected at three distinct moments at 600, 900, and 1200 ms (expected signals, ESs). Here the RS occurs finally after the longest possible delay of 1500 ms. After the RS, the requested movement was executed (reaction and movement time, RT-MT). See [8,84] and the text for the detailed experimental protocol. (B) Dot displays of the spike sequences (duration: 2 s, sampling resolution 1 ms) of three neurons simultaneously recorded from the primary motor cortex (MI). The spike synchronous spike events across the three neurons to ESs is no longer apparent. However, the state-space log-linear model reveals a gradual increase in the strength of the triple-wise interaction until the end of the PP.
To strengthen the findings that a triple-wise interaction term increases during the PP, we test for the presence of a triple-wise spike correlation in the PP using the Bayes factor (marginal likelihood ratio, Eq. 12). We compute the BF for the opposing models using a full model (r~3) of N~3 neurons for the hidden log-linear parameters. In one hypothetical model, M 1 , we assume that the triple-wise interaction term is positive, S 1 : h t 123 w0, in the other model, M 2 we assume that the triple-wise interaction term is smaller than or equal to zero, S 2 : h t 123 ƒ0. A large positive value of the BF indicates that the spike data in that period supports model M 1 as opposed to model M 2 , i.e., the existence of excess synchronous spike events among the three neurons in comparison with the chance rate given by spike rates and pairwise correlations. A large negative value of the BF shows support of model M 2 as opposed to model M 1 , i.e., a paucity of synchronous spikes for the three neurons. Figure 7A shows the BF computed bin-by-bin (Eq. 13) in bit units. In the bottom panel of Figure 7A, we indicate the behavioral periods for which we test for the presence of triple-wise spike correlation (the PP and RT to MT; in later analysis we divide the PP into early and late stages). The evidence for model M 1 , as opposed to M 2 , in the behaviorally relevant time periods is obtained by summing the log of the bin-by-bin BF in that period (cf. Eqs. 12 and 13). The BF in the PP is found to be 18.08 bit, which is interpreted as 'very strong' evidence for presence of a positive triple-wise spike correlation by the classical guideline [76].
Next, we test the observed value of the BF in the PP (18.08 bit) by comparing it with surrogate BFs. These are computed from null data that contain no triple-wise spike interactions, while keeping the timevarying structure of the pairwise correlations and the individual spike rates the same as those observed in the original spike train data. The construction of the null data follows the same procedure as in the simulation study shown before: We apply a pairwise state-space log-linear model (r~2) to the spike data and then generate 1000 surrogate samples (each with n~36 trials) of N~3 parallel spike sequences using the fitted pairwise state-space log-linear model. Figure 7B (left panel) displays the surrogate distribution along with the observed BF in the PP. The observed positive BF falls out of the 95% confidence interval, suggesting the presence of a positive triplewise interaction as an underlying model for the spike train data in this period. In the right panels of Figure 7B, we display the results of the same analysis, but using spike train data binned by a larger (D~5 ms, top panel) and a smaller (D~2 ms, bottom panel) binwidth. We obtained almost the same results as in the analysis with a bin-width of 3 ms as for the analysis with the larger bin-width (5 ms). However, we could not reach the same results for the smaller bin width (2 ms) because the small count of synchronous spike events made it impossible to reject the null hypothesis.
If the higher-order dependency among the three neurons is related to motion preparation, the evidence for the triple-wise spike correlation should be stronger in the late stage of the PP than before. To test this idea, we divided the PP into earlier and later stages of the PP (each with a 750 ms duration, see the bottom of Figure 7A), and investigate whether the three MI neurons exhibit a triple-wise spike correlation in these periods. In addition, we select a duration of 300 ms after the onset of the RS, during which the animal starts to initiate the motor action (reaction and movement time, RT-MT). The upper panels of Figure 7C display the results obtained with bins of width D~3 ms. The weights of evidence for the triple-wise spike correlation models are 7.16, 11.2, and 0.98 bit in the early PP, late PP, and RT-MT periods respectively: i.e., 'very strong' evidence is found at the late stage of the PP. The null hypothesis is rejected in this late period of the PP whereas the same null hypothesis is not rejected in the earlier period of PP or the period after the RS. In the late PP, the significantly large positive BF indicates the existence of a positive triple-wise spike correlation. The bottom panels of Figure 7C display the same analyses, but with a larger (D~5 ms) bin-width. We obtain the same results as in the analysis with a bin-width of 3 ms. When we analyze the data using a smaller bin-width (D~2 ms), the BFs are not significant in any of the periods (not shown here) because of the small samples of synchronous spikes. In addition, as in the simulation study in the previous subsection, we tested whether the observed synchronous spiking activities are explained merely by simultaneous increases in the pairwise interaction terms of the second order log-linear model. We did not find any evidence for such simultaneous increases in the pairwise interactions from the data for all three periods (bin-width D~3 ms) (not shown here). Thus only by the application of the higher-order analysis we were able to detect the task-dependent changes in the joint interactions of all three neurons.
The critical assumptions made in this study are independence and identical sampling across the trials (across-trial stationarity). However, in the spike sequences of neuron 2 (neuron id 3) shown in Figure 6B, we notice an increase in the firing rates across trials. Higher firing rates yield a higher chance rate for synchronous spiking events. Thus, underestimation of the spike rate caused by averaging across the entire trials might induce spurious estimation of higher-order spike interaction. Indeed, synchronous spike events of the three neurons are mostly observed in later trials (e.g., trials 19-36, see blue dots in Figure 6B) when the firing rates were higher than in the earlier trials. Hence, we repeat the same analysis by using only the latter half of the trials (trials [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36]. Figure 7D displays the results. We still observe a significantly positive BF in the late PP. We note that the estimated dynamics of the log-linear parameters using only the latter half of the trials are not changed much from those using all trials, i.e., Figure 6D. In the analysis using the first half of the trials (1-19, not shown), the BF is not significant during the late PP. Kilavik et al. [19] reported that ''during practice, the temporal structure of synchrony was shaped, with synchrony becoming stronger and more localized in time during late experimental sessions''. We observe similar effects in our triple-wise spike correlation analysis.  Table 1. The panels depict the MAP estimates of the log-linear interaction parameters (solid lines; from top to bottom, the first and second order log-linear parameters and a triple-wise spike interaction, h t 123 ). The gray bands indicate the 95% credible interval. In this analysis, we used an identity matrix as an autoregressive parameter, i.e. F~I, in the state model. The covariance matrix of the initial parameter was fixed to S~0:05I. (E) The smoothed estimate of a triple-wise spike interaction, h t 123 , is computed from binary data constructed using a bin-width of 2 ms (Top) and a bin-width of 5 ms (Bottom). The top of each panel shows the timing of the synchronous spike events of all three neurons. doi:10.1371/journal.pcbi.1002385.g006

Discussion
In this study, we introduced a novel method for estimating dynamic spike interactions in multiple parallel spike sequences by means of a state-space analysis (see Methods for details). By applying this method to nonstationary spike train data using the pairwise log-linear model, we can extend the stationary analysis of the spike train data by the Ising/spin-glass model to within-trial nonstationary analysis (Figure 3). In addition, our approach is not limited to a pairwise analysis, but can perform analyses of timevarying higher-order spike interactions (Figure 4). It has been discussed whether higher-order spike correlations are important to characterize neuronal population spiking activities, assuming stationarity in the spike data [54][55][56][57][58][59]. Based on the state-space model optimized by our algorithm, we developed two methods to validate and test its latent spike interaction parameters, in particular the higher-order interaction parameters, which may dynamically change within an experimental trial. In the first method, we selected the proper order for the spike interactions incorporated in the model under the model selection framework using the approximate formula of the AIC for this state-space model (In Methods, 'Selection of state-space model by information criteria'). This method selects the model that best fits the data overall across the entire observation period. The selected model can then be used to visualize the dynamic spike interactions or for a performance comparison with other statistical models of neuronal spike data. However, more importantly, the detailed structure of the transient higher-order spike interactions needs to be tested locally in time, particularly in conjunction with the behavioral paradigm. To meet this goal, we combined the Bayesian model comparison method (the Bayes factor) with a surrogate method (In Methods, 'Bayesian model comparison method for detecting spike correlation'). The method allows us to test for the presence of higher-order spike correlations and examine its relations to experimentally relevant events. We demonstrated the utility of the method using neural spike train data simultaneously recorded from primary motor cortex of an awake monkey. The result is consistent with, and further extended the findings in the previous report [8]: We detected an increase in triple-wise spike interaction among three neurons in the motor The bin-by-bin BF is computed from Eq. 13, using the state-space loglinear model fitted to the spike data (a total of 36 trials, 2 s binned using 3 ms bin-width; cf. Figure 6C). The BF computes evidence of a positive triple-wise spike interaction, M 1 : h t 123 w0, as opposed to a zero or negative triple-wise interaction, M 2 : h t 123 ƒ0. The evidence for model M 1 as opposed to M 2 in a behaviorally relevant time period is obtained by summing the log of the bin-by-bin BF in that period (cf. Eq. 12). (Bottom) Two behavioral periods (preparatory period, PP, and reaction and movement time, RT-MT) tested for presence of a triple-wise spike correlation. In addition, to examine the evolution of the triple-wise spike correlation in the PP, the PP is divided into early and late stages at the middle of the PP. (B) (Left) The observed BF for an entire PP (marked by a red line and triangle) computed using Eq. 12 (bin-width: 3 ms). We then test the 'observed BF' using a distribution function (solid line) of null BFs derived from surrogate data sets generated from a fitted model containing only up to pairwise interaction terms (r~2). The gray area indicates the 95% confidence interval of the distribution. Vertical cortex during the preparatory period in a delayed motor task, which was also tightly locked to the expected signals. Although the analysis was done for a limited number of neurons, smaller than the expected size of an assembly, it demonstrates that the nonstationary analysis of the higher-order activities is useful to reveal cooperative activities of the neurons that are organized in relation to behavioral demand. Of course, further analysis is required to strengthen the findings made above including a metaanalysis of many different sets of multiple neurons recorded under the same conditions.
In this study, we adopted the log-linear model to describe the higher-order correlations among the spiking activities of neurons. There are, however, other definitions for 'higher-order spike correlation'. An important alternative concept is the definition based on cumulants. Using the cumulants of an observed count distribution from a spike train pooled across neurons, Staude et al. developed an iterative test method that can detect the existence of a high amplitude in the jump size distribution of the assumed compound Poisson point process (CPP) model for the pooled spike train [34,35]. This method can detect an assembly from a few occurrences of synchronous spike events to which many neurons belong to, typically by using the lower-order cumulants of the observed spike counts. In contrast, the information geometry measure for the higher-order spike correlation used in this study aims to represent the correlated state that cannot be explained by lower-order interactions. Consequently, the information geometry measure extracts the relative strength of the higher-order dependency to the lower-order correlated state. Therefore, the presence of positive higher-order spike correlations does not necessarily indicate that many neurons spike synchronously whenever they spike because such activities can be induced by positive pairwise spike correlations alone [65,85,86] (see also Figure 8A and B in the Methods section). In contrast, the cumulant-based correlation method by Staude et al. [34] infers the presence of 'higher-order correlation' for such data by determining the presence of high amplitudes in the jump size distribution of the assumed CPP model. Yet another important tool for analyzing higher-order dependency among multiple neurons is the copula function, a standardized cumulative distribution function used to model the dependence structure of multiple random variables (see [87][88][89] for an analysis of neurophysiological data using the copula, including an analysis for higher-order dependency [89]). In summary, it should be remembered that the analysis method used for the higher-order dependency among neuronal spikes inherits its goal from the assumed model for spike generation as well as a parametric measure defined for the 'higher-order' spike correlation [34,62].
Although we face a high-dimensional optimization problem in our settings, we are able to successfully obtain MAP estimates of the underlying parameters because of the simplicity of the formulation of the state-space model: The use of the log-concave exponential family distributions [50,90] in both the state and observation models guarantees that the MAP estimates can be obtained using a convex optimization program. At each bin, the method numerically solves a nonlinear filter equation to obtain the mode of the posterior state density (the MAP estimates, see Eqs. 28,29,and 30 in Methods). With only a few (3-8) Newton-Raphson iterations, the solution reaches a plateau (the increments of all the elements of the updated state space vector are smaller than 10 {5 ). The entire optimization procedure can be performed in a reasonable amount of time: On a contemporary standard laptop computer it takes no longer than 30 seconds to obtain smooth estimates of a full log-linear model for N~3 neurons (T~500 bins, Figure 4), which includes 100 EM iterations. The method is even faster when approximating the posterior mode using the update formula Eq. 30 without any iterations, using the one-step prediction mean as an initial value. This fast approximation method suggested in [69] could even be utilized in a realtime, on-line application of our filter (the filtering method applied to a single trial, n~1, using predetermined hyper-parameters) at the cost of estimation accuracy.
The pairwise analysis can be applied up to about N~12 neurons simultaneously to derive time-dependent pair interactions. However, the current version of the algorithm does not scale to a larger number of neurons because the number of spike patterns that need to be considered suffers from a combinatorial explosion. The major difficulty arises from the coordinate transformation from the h-coordinates to the g-coordinates that appear in the non-linear filter equation (Eq. 31 in Methods). The coordinate transformation is required in this equation to calculate the innovation signal, i.e., the difference between the observed synchrony rates, y t , and the expected synchrony rates (gcoordinates) based on the model. We numerically derived the exact g-coordinates by marginalizing the 2 N dimensional joint probability mass function computed from the h-coordinates. Thereby, a full knowledge of the probability mass function is required even if the model considers only the lower-order interactions. Because this is a common problem in the learning of artificial neural networks [46,47,91], sampling algorithms such as the Markov chain Monte Carlo method have been developed to approximate the expectation parameters, h, without having to compute the partition function [92]. The inclusion of such methods allows us to analyze the time-varying low order spike interactions from a larger number of parallel spike sequences. Recent progress [93], e.g., in the mean field approach and/or the minimum probability flow learning algorithm for an Ising model, may allow us to further increase the number of neurons that can be treated in this nonstationary pairwise analysis. Nonetheless, the method presented here, which aims at a detailed analysis of the dynamics in higher-order spike interactions, may not easily scale to massively parallel spike sequences that can be analyzed by other methods such as those based on the statistics pooled across neurons. Thus, we consider it to be important to combine the detailed analysis method proposed in this study with other state-ofthe-art analysis techniques in practical applications. For example, test methods based on population measures such as the Unitary Event method and cumulant-based inference method [34,35] allow us to detect the existence of statistically dependent neurons in massively parallel spike sequences. If the null-hypothesis of independence among those neurons is not rejected in these methods, we can exclude those neurons from any further detailed analysis of their dynamics using the methods proposed in this study.
Several critical assumptions made in the current framework need to be addressed. First, it was assumed in constructing the likelihood (Eq. 7) that no spike history effect exists in the generation of a population spike pattern. Second, we assumed the use of identically and independently distributed samples across trials when constructing the likelihood (Eq. 7). The first assumption may appear to be strong constraint given the fact that individual neurons exhibit non-Poisson spiking activities [94]. However, as in the case of the estimation of the firing rate of a single neuron, the pooled spike train across the (independent) trials is assumed to obey a Poisson point process because of the general limit theorem for point processes [30,31,95,96]. This is because most of the spikes in the pooled data come from independent different trials. They are thus nearly statistically independent from each other, even if the individual processes are non-Poisson. Similarly, in our analysis, we used statistics from a pooled binary spike train, assuming independence across trials: The occurrences of joint spikes in the binary data pooled across trials are mostly independent of each other across bins. Because these joint spike occurrences are sparse (i.e., they rarely happen closely to each other in the same trial), it is even more feasible to assume their statistical independence across bins. Third, however, while pooling independent and identical trials (the second assumption) may validate the first assumption of the independence of the samples across bins, that assumption of independently and identically distributed samples across trials has itself been challenged [97,98] and is known to be violated in some cases, e.g., by drifting attention, ongoing brain activity, adaptation, etc. It is possible that the trial-by-trial jitter/variation in the spike data causes spurious higher-order spike correlation. Thus, as discussed in the section on the application of our methods to real neuronal data, it is important to examine the stationarity of the spike train data across trials. Note that, not only the firing rates, but also the spike synchrony can be shaped on a longer time scale by repeatedly practicing a task [19]. In fact, the current analysis method can be used to examine the long-term evolution of pairwise and higherorder spike interactions across trial sessions by replacing the role of a bin with a trial, assuming within-trial stationarity. It will be a challenge to construct a state-space log-linear model that additionally applies a smoothing method across trials (see [98] for such a method for a point process model).
The present method is left with one free parameter, namely the bin-width D. The bin-width determines the permissible temporal precision of synchronous spike events. Very large bin-widths result in binary data that are highly synchronized across sequences, while very small bin-widths result in asynchronous multiple spike sequences. In the latter case, we might overlook the existing dependency between multiple neural spike sequences due to disjunct binning [99] (but see [57,100] that aim to overcome such a problem by modeling the spike interactions across different consecutive time bins). Within our proposed modeling framework, which focuses on instantaneous higher-order spike correlations, it is important to catch the innate temporal precision of the neuronal population under investigation using the appropriate bin-width. Thus, the choice may be guided by the biophysical properties of the neurons. However, it may be of advantage to derive the binwidth in a data-driven manner. For example, in the context of an encoding problem, the proper bin-width can be chosen based on the goodness-of-fit test for single neuron spiking activities [101], conditional on the spiking activities of the other neurons [40]. For questions about the relation of coincident spiking to stimulus/ behavior, the bin-width may be selected based, for example, on the predictive ability of an external signal. For this goal, it is important to search the optimal bin-width using elaborate methods such as those developed in the context of the Unitary Event analysis method [99] (see [25] for a review of related methods).
A substantial number of studies have demonstrated that stimulus and behavioral signals can be decoded simply based on the firing rates of individual neurons. At the same time, it has been discussed whether spike correlations, particularly higher-order spike correlations, are necessary to characterize neuronal population spiking activities [54][55][56][57][58] or to encode or decode information related to stimuli [53,60,102]. At this point in time, a smaller number of dedicated experiments have supported the conceptual framework of information processing using neuronal assemblies formed by neurons momentarily engaged in coordinated activities, as expressed by temporally precise spike correlations (see [6,7,9,10] for reviews of these experiments). Nevertheless, it is possible that the current perspective on this subject has been partly formed by a lack of proper analysis approaches for simultaneously tracing time-varying individual pairwise spike interactions, and/or their higher-order interactions. Indeed, we demonstrated by the time-resolved higher-order analysis that three cortical neurons coordinated their spiking activities in accordance with behaviorally relevant points in time. Thus our suggested analysis methods are expected to be useful to reveal the dynamics of assembly activities and their neuronal composition, as well as for testing their behavioral relevance. We hope that these methods help shed more light on the cooperative mechanisms of neurons underlying information processing.

Mathematical properties of log-linear model
In this subsection, we review the known mathematical properties of a log-linear model for binary random variables. These properties are used in constructing recursive filtering/ smoothing formulas in the next section. Using the multi-index, I[fV 1 , . . . ,V N g (see the Results subsection 'Log-linear model of multiple neural spike sequences'), the probability mass function (Eq. 1), p(x), where x~x 1 ,x 2 , . . . ,x N ð Þ and x i~1 or 0 (i~1, . . . ,N), and the expectation parameters (Eq. 2) are compactly written as and where f I x ð Þ is a feature function, here representing an interaction among the neurons indicated by the multi-index, I (Eq. 3).
The hand g-coordinates are dually flat coordinates in the exponential family probability space [44,50], and the coordinate transformation from one to the other is given by the Legendre transformation [40,50]. From Eq. 14, the log normalization function, y(h), is written as The first derivative of the log normalization function, y(h), with respect to h I (I[fV 1 , . . . ,V N g), provides the expectation parameter, g I~E ½f I (X): Let Q(h) be the negative entropy of the distribution: Eqs. 17 and 18 complete the Legendre transformation from hcoordinates to g-coordinates [44,50]. The Legendre transformation transfers the functional relationship of h and y(h) to the equivalent relation in the dual coordinates, g and Q(g). The inverse transformation is given by Eq. 18 and LQ(g) Lg I~h I .
Using the log normalization function, we can obtain the multivariate cumulants of p(x) with respect to the random variables, f I (X). The cumulant generating function of the exponential family distribution is given as K u ð Þ~y(hzu) {y(h). Let us compactly write the partial derivative with respect to h I (i.e., L Lh I ) as L I . Then, the first order cumulant is given as L Lu K u ð Þ u~0~L I y(h)~g I , as shown in Eq. 17. In general, the cumulants of the exponential family distribution are given by the derivatives of the log normalization function. Thus, the second derivative of y(h) yields the second-order cumulant, g IJ :E½(f I (X){g I )(f J (X){g J )~g I|J {g I g J (by the cup, |, we mean the multi-index representation of an union of the elements of the two multi-indices, e.g., if I~12 and J~234, then I|J~1234): for I,J[fV 1 , . . . ,V N g. g IJ is known as the Fisher metric with respect to the natural parameters. Eqs. 17 and 19 are important relations used in this study to construct a non-linear filtering equation for a dynamic estimate of the natural parameters because we approximate the log-linear model (Eq. 14) with a precision of up to a (log) quadratic function (cf. Eqs. 28 and 29). Similarly, the higher-order derivatives yield higher-order multivariate cumulants. For example, the third-order derivative yields the third order cumulant, L I L J L K y(h)~C ? IJ,K , where C ?
The pseudo distance between two different distributions, q(x) and p(x) is defined using the Kullback-Leibler (KL) divergence We represent distribution q(x) by using g-coordinates as g q ð Þ, and p(x) by using h-coordinates as h p ð Þ. Here, we used h p ð Þ for the hparameters of p(x) (and g q ð Þ for g-parameters of q(x)) in order to differentiate it from the representation of distribution q(x) in the hcoordinates (and the representation of p(x) in the g-coordinates). Then, the KL-divergence between the two distributions, q(x) and p(x), is computed as [44,50] D½q,p~y(h p ð Þ)zQ(g q ð Þ){h p ð Þg q ð Þ: ð21Þ Optimized estimation of dynamic spike interactions We develop a non-linear recursive Bayesian filtering/smoothing algorithm and its optimization method in order to trace dynamically changing spike interactions from parallel spike sequences. To reach this goal, we use the expectation-maximization (EM) algorithm [68,73,103,104], which is known to efficiently combine the construction of the posterior density of a state (the natural parameters) and the optimization of the hyper-parameters. This method maximizes the lower bound of the log marginal likelihood, Eq. 10. Using Jensen's inequality and nominal hyperparameters, w~F,Q,m,S ½ , the lower bound of the log marginal likelihood with hyper-parameters w Ã~FÃ ,Q Ã ,m Ã ,S Ã ½ is given by Here, H represents a negative entropy. The maximization of the lower bound with respect to w Ã is equivalent to maximizing the expected complete data log-likelihood in Eq. 22, known as the Qfunction: Q w Ã jw ð Þ~E log p y 1:T ,h 1:T jw Ã ð Þ j y 1: The expectation in the above equation is read as E½ jy 1:T ,w. We maximize the Q-function by alternating the expectation (E) and maximization (M) steps. In the E-step, we obtain the expected values with respect to h t in Eq. 23 using a fixed w. In the M-step, we obtain the hyper-parameter, w Ã , that maximizes Eq. 23. The resulting w Ã is then used in the next E-step. The details of each step are now given as follows. E-step: Bayesian recursive filter/smoother. The E-step is composed of filtering and smoothing algorithms conducted by forward and backward recursions, respectively. The forward algorithm sequentially constructs a posterior density of the state at time t given the spike data up to and including time t, whereas the backward algorithm constructs a posterior density at time t given the entire data. The posterior density allows us to compute the maximum a posteriori (MAP) estimate or Bayes estimator and provides uncertainty for the estimate. In the following, h tjs and W tjs denote the conditional mean, E h t jy 1:s ,w ½ , and covariance, E½(h t {h tjs )(h t {h tjs ) 0 jy 1:s ,w. The filter mean and covariance are denoted as h tjt and W tjt , respectively. The mean and covariance of a smooth posterior density are denoted as h tjT and W tjT , respectively.
We first compute the one-step prediction density, p h t jy 1:t{1 ,w ð Þ . This is the conditional density of the state at time t given the observation of parallel spike sequences up to time t{1. The one-step prediction density is written using the Chapman-Kolmogorov equation as [37,69,105] p h t jy 1:t{1 ,w ð Þ Here the transition probability, p h t jh t{1 ,w ð Þ , is a multivariate normal distribution with mean Fh t{1 and covariance Q, as defined in the state equation, Eq. 8. For an initial prior, p h 1 ð Þ, we use a normal distribution with mean h 1j0~m and covariance W 1j0~S . The other distribution, p h t{1 jy 1:t{1 ,w ð Þ , in Eq. 24 is the filter density at time t{1. In the next paragraph, the filter density will be obtained by approximating it with a normal distribution whose mean and covariance are denoted as h t{1jt{1 and W t{1jt{1 . Under this condition, the one-step prediction density (Eq. 24) again becomes a normal distribution whose mean, h tjt{1 , and covariance, W tjt{1 , are given as [37,68] h tjt{1~F h t{1jt{1 , ð25Þ The filter density, p h t jy 1:t ,w ð Þ , is the conditional distribution of the state given the observation of parallel spike sequences up to time t. Using the likelihood function and one-step prediction density, the filter density is given by Bayes' theorem as p h t jy 1:t ,w ð Þ p y t jh t ,y 1:t{1 ð Þ p h t jy 1:t{1 ,w ð Þ p y t jy 1:t{1 ,w ð Þ This posterior density is a complicated function with respect to the natural parameter, h t . Here, we apply a Gaussian approximation to the posterior density using the Laplace method [37,104,106,107]: The filter mean, h tjt , is identified with a mode of the posterior density as h tjt~a rgmax q log p hjy 1:t ,w ð Þ , and the filter covariance, W tjt , is determined from the Hessian of the log posterior at the mode as W tjt~{ ½++log p hjy 1:t ,w ð Þ j h~h tjt {1 .
The approximate posterior mode is obtained using the iterative procedure of a gradient ascent method or the Newton-Raphson method using the gradient and Hessian matrix. The gradient and Hessian of the log of the posterior density at h is calculated as +f:+log p hjy 1:t ,w The gradient and Hessian, Eqs.28 and 29, are evaluated using an old value, h old . Here, E is a learning coefficient that was introduced in the context of 'natural' gradient search algorithm [91]. For the small population size analyzed in this study, we used E~1. However, it is recommended that a smaller positive value be selected for the analysis of a larger system to avoid numerical instability. Because both the likelihood function and prior density are logarithmically concave functions, the posterior is also logconcave. Thus, this optimization problem is convex, which guarantees a unique solution for the filter estimate [90]. The optimized natural parameter is selected as the filter mean, h tjt . The filter covariance is approximated as W tjt~{ H {1 by using h tjt . It is also possible to use a simple gradient ascent method to obtain the mode, and then compute the Hessian matrix at the mode. The above method is equivalent to solving the following nonlinear recursive filter formulas: Eq. 31 was obtained from +log p h t jy 1:t ,w ð Þ0. Eqs. 31 and 32 are recursively computed in combination with the one-step prediction equations, Eqs. 25 and 26, for t~1, . . . ,T. Because h tjt and g tjt are dual representations of the same probability distribution, Eq. 31 is a nonlinear equation and needs to be solved as suggested above (Eqs. 28,29,30). As pointed out for a point process adaptive filter [37,69], the residual, y t {g tjt , in Eq. 31 acts similarly to an innovation vector of a standard Kalman filter. The same error signal between the observed synchrony rates and expected synchrony rates is utilized in training the Boltzmann machine [46,47,91]. The innovation term corrects the one-step prediction mean, h tjt{1 , i.e., an expected state by the prior distribution. The degree of the correction is determined by the number of repeated trials, n, and uncertainty of the predicted state, W tjt{1 . The hyper-parameters of the prior density significantly affect the latter gain (see Eq. 26), and thereby the smoothness of the estimated processes. The filter covariance equation, Eq. 32, describes the reduction of the prediction uncertainty, W tjt{1 , by observing the parallel spike sequences at time t, with the amount determined by the number of repeated trials, n, and the Fisher information, G tjt . Please see Figure 1 for a geometric view of the recursive Bayesian filter.
Finally, we compute the smooth posterior density using the backward recursive algorithm [68,105,107], p h t jy 1:T ,w ð Þp h t jy 1:t ,w ð Þ ð p h tz1 jy 1:T ,w ð Þ p h tz1 jh t ,w ð Þ p h tz1 jy 1:t ,w ð Þ dh tz1 : ð33Þ Because the density functions in the recursive formula were approximated as a normal distribution, we follow the fixed-interval smoothing algorithm [37,68,105] established for a Gaussian state and observation equation. Starting from h TjT and W TjT , which are obtained from the filtering algorithm, we obtain the smoothed mean and covariance, with for t~T{1,T{2, . . . 1. The lag-one covariance smoother, W t{1,tjT , which appears in the Q-function, is obtained using the method of De Jong and Mackinnon [108]: M-step: Optimization of hyper-parameters. At the Mstep, we optimize hyper-parameter w given the posterior density under the principle of maximizing the Q-function. From L LQ Ã Q w Ã jw ð Þ~0, the update rule of the covariance matrix, Q, is obtained as Here, h tjT , W tjT , and W t{1,tjT are the smoother mean and covariance, and the lag-one covariance matrix given by Eqs. 34, 35, and 37, respectively. Similarly, from L LF Ã Q w Ã jw ð Þ~0, the auto-regressive parameter, F, is updated according to The mean of the initial distribution is updated with m Ã~h 1jT from L Lm Ã Q w Ã jw ð Þ~0. The covariance matrix, S, is not updated; instead, we use a fixed matrix, S. Alternatively, the covariance matrix of the initial distribution can be updated according to S Ã~W 1jT z(h 1jT {m)(h 1jT {m)' from L LS Ã Q w Ã jw ð Þ~0, while the mean, m, is fixed. Both methods work well with appropriate choices for the fixed parameters. In this study, we updated the mean vector, m, of the initial normal distribution, and used a fixed diagonal matrix for its covariance matrix, S. It was also suggested to use a stationary mean and covariance of an unconstrained state process (Eq. 8) as parameters of the initial prior distribution [68,109]. The equilibrium condition from Eq. 8 yields m~Fm and S~F X F 0 zQ. The solutions are obtained as m Ã~0 and vec(g Ã )~½I{F6F {1 vec(Q) (p. 121 and p. 426 in [109], p. 112 in [110]). In the latter, 6 denotes the Kronecker product (tensor product) and the vec operator creates a single column vector from a matrix by stacking its column vectors. However, the stationary condition is not always satisfied. As we found that the fitted F sometimes violates the stationary assumption, we did not adopt this approach in the current study.

Selection of state-space model by information criteria
The method developed in the previous subsection is applicable to a full log-linear model, as well as a model that considers an arbitrary order of interactions. In order to select the most predictive model among the hierarchical models in accordance with the observed spike data, we select the state-space model that minimizes the Akaike information criterion (AIC) for a model with latent variables [73,105]: Here, l w ð Þ is the log marginal likelihood (Eq. 10) and k is the number of free hyper-parameters in the prior distribution. For the rth order model, the number of natural parameters is given by d~X r k~1 N k . The number of free parameters in the prior distribution is computed as k~d 2 zd(dz1)=2zd, where each term corresponds to the number of free parameters in F, Q, and m. Note that the AIC applied to the state-space model is sometimes referred to as the Akaike Bayesian information criterion (ABIC) [73]. In the following, we derive the approximation method to evaluate the AIC for the state-space log-linear model.
The log marginal likelihood, Eq. 10, can be written as l(w)~X T t~1 log p(y t jy 1:t{1 ,w) X T t~1 log ð p(y t jh t )p(h t jy 1:t{1 ,w)dh t : We make a log quadratic approximation to evaluate the integral.
To accomplish this, we denote ð p(y t jh t )p(h t jy 1:t{1 ,w)dh t~1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi with q(h t )~n y The Laplace approximation of the integral in Eq. 42 is given as [80] ð exp q h t We confirmed that the log-quadratic approximation provided a better estimate of the marginal likelihood than the first order approximation used in [77] by comparing them with a Monte Carlo approximation of the integral in Eq. 42. We select the statespace log-linear model that minimizes the AIC (Eq. 40), where the log marginal likelihood is approximated using Eq. 45. For comparison with the AIC, we compute two other information criteria that employ different forms of the penalization term. The Bayesian information criterion [79,80] (also known as Schwartz's criterion) are obtained by replacing the penalization term of Eq. 40, k, with k log n: BIC~{2 l(w)zk log n: Shimodaira's predictive divergence for indirect observation models (PDIO) [81] is given as Here, we redefine w as a one-dimensional vector of free hyperparameters, while M(w) denotes the one-step operator of EM iteration. To obtain the Jacobian matrix for the EM operator, we follow the algorithm described in Meng and Rubin [111]. In this method, the Jacobian matrix was approximated using a numerical differentiation of the EM operator. By perturbing one hyperparameter and then computing a one-step EM iteration, numerical differentiations of the hyper-parameters with respect to the perturbed hyper-parameter were obtained. An entire Jacobian matrix was approximated by repeating the process while changing the hyper-parameter to be perturbed.

Bayesian model comparison method for detecting spike correlation
In this subsection, we formulate a method for detecting the hidden structure of spike interaction by means of a Bayesian model comparison based on the Bayes factor (BF) [74][75][76]. The BF is a ratio of likelihoods for the observed data, based on two different assumptions about the hidden states (model M 1 and M 2 ). Here we reiterate the definition of the BF for model M 1 as opposed to model M 2 used in this paper (cf. Eq. 12): B 12 (y a:b )~p y a:b jM 1 ,w ð Þ p y a:b jM 2 ,w ð Þ : The BF becomes larger than 1 if the data, y a:b in a time period a,b ½ , supports model M 1 as opposed to model M 2 , and becomes smaller than 1 if the data supports model M 2 as opposed to model M 1 . The BF can be computed from the one-step prediction and triple-wise correlation as evaluated by Test 2. The existing positive triple-wise interaction present in period III is well detected by the large BF, and the very low BFs in period I and II correctly detect the absence of triple-wise correlation.
In Figure 8D, we examine how much each individual spike pattern contributes to the calculation of the BFs discussed above. The weight of evidence (log of the BF) in an observation period ½a,b is obtained by summing the pieces of evidence computed locally at time t using Eq. 52, as log 2 B 12 (y a:b )P b t~a log 2 B 12 (y t jy a:t{1 ). Here we elucidate the contributions of each individual spike pattern to the BF by sorting the bin-by-bin evidence, log 2 B 12 (y t jy a:t{1 ), with each spike pattern. Figure 8D displays the average BFs for each spike pattern in the three periods in bit unit. The degree of the contribution by a specific spike pattern to the evidence for the tested spike interaction model varies depending on the context in which the spike pattern is observed. For example, the magnitudes of the BFs for a triple-wise correlation (Test 2) found by observing a spike triplet (111) vary greatly across the three periods (see dark gray bars at 111 in periods I, II, and III). In period III, where pair-synchronous spikes are rarely observed because of the existence of negative pairwise interaction terms, the observation of triplet spikes provides substantial evidence of a triple-wise correlation. In contrast, in period II, where pair-synchronous spikes frequently occur because of the existence of positive pairwise interactions, the chance triplet spikes do not provide substantial evidence of a triple-wise correlation. These results come from the fact that an unexpected synchronous spike pattern that significantly alters and updates the filter density from its one-step prediction density gives rise to a large absolute BF value (Eq. 52). By the same logic, the BFs for the 000 and 100 patterns (and 010 and 001, not shown) are small because they do not substantially change the posterior densities. However, they should not be neglected because these patterns are abundant as a result of the low firing rates. For example, in period I in Figure 8C, the accumulation of the small weight of evidence from the abundant spike patterns (such as 000 and 100) offsets the large weight of evidence induced by a few chance coincidences such as triplet spikes (111), thus producing virtually zero weight of evidence for this period.