Single-trial estimation of stimulus and spike-history effects on time-varying ensemble spiking activity of multiple neurons: a simulation study

Neurons in cortical circuits exhibit coordinated spiking activity, and can produce correlated synchronous spikes during behavior and cognition. We recently developed a method for estimating the dynamics of correlated ensemble activity by combining a model of simultaneous neuronal interactions (e.g., a spin-glass model) with a state-space method (Shimazaki et al. 2012 PLoS Comput Biol 8 e1002385). This method allows us to estimate stimulus-evoked dynamics of neuronal interactions which is reproducible in repeated trials under identical experimental conditions. However, the method may not be suitable for detecting stimulus responses if the neuronal dynamics exhibits significant variability across trials. In addition, the previous model does not include effects of past spiking activity of the neurons on the current state of ensemble activity. In this study, we develop a parametric method for simultaneously estimating the stimulus and spike-history effects on the ensemble activity from single-trial data even if the neurons exhibit dynamics that is largely unrelated to these effects. For this goal, we model ensemble neuronal activity as a latent process and include the stimulus and spike-history effects as exogenous inputs to the latent process. We develop an expectation-maximization algorithm that simultaneously achieves estimation of the latent process, stimulus responses, and spike-history effects. The proposed method is useful to analyze an interaction of internal cortical states and sensory evoked activity.


Introduction
Neurons in the brain make synaptic contacts to each other and form specific signaling networks. A typical cortical neuron receives synaptic inputs from 3000 − 10000 other neurons, and makes synaptic contacts to several thousands of other neurons. They send and receive signals using pulsed electrical discharges known as action potentials, or spikes. Therefore, individual neurons in a circuit can be activated in a coordinated manner when relevant information is processed. In particular, nearly simultaneous spiking activity of multiple neurons (synchronous spikes) occurs dynamically in relation to a stimulus presented to an animal, the animal's behavior, and the internal state of the brain (attention and expectation) [1][2][3][4][5].
Recently, it was reported that a model of synchronous spiking activity that accounts for spike rates of individual neurons and interactions between pairs of neurons can explain ∼ 90% of the synchronous spiking activity of a small subset (∼ 10) of retinal ganglion cells [6,7] and cortical neurons [8] in vitro. This model is known as a maximum entropy model or an Ising/spin-glass model in statistical physics. However, since the model assumes stationary data, it is not directly applicable to non-stationary data recorded from awake behaving animals. In these data sets, spike-rates of individual neurons and even interactions among them may vary across time.
In order to analyze time-dependent synchronous activity of neurons, we recently developed a method for estimating the dynamics of correlations between neurons by combining the model of neuronal interactions (e.g., the Ising/spin-glass model) with a state-space method [9,10]. In classical neurophysiological experiments, neuronal activity is repeatedly recorded under identical experimental conditions in order to obtain reproducible features in the spiking activity across the 'trials'. Typically, neurophysiologists estimate average time-varying firing rates of individual neurons in response to a stimulus from the repeated trials [11,12]. In the same fashion, the statespace method in [9,10] aims to estimate the dynamics of the neuronal interactions, including higher-order interactions, that occurs repeatedly upon the onset of externally triggered events. When this method is applied to three neurons recorded simultaneously from the primary motor cortex of a monkey engaged in a delayed motor task (data from [2]), it was revealed that these neurons dynamically organized into a group characterized by the presence of a higher-order (triple-wise) interaction, depending on the behavioral demands to the monkey [12].
However, neurophysiological studies in the past decades revealed that spiking activity of individual neurons is subject to large variability across trials due to structured ongoing activity of the networks that arises internally to the brain [13][14][15]. In these conditions, the method developed in [12] would not efficiently detect the stimulus responses because a signal-to-noise ratio may be small even in the trial-averaged activity. Although statistical methods for detecting responses of individual neurons from single-trial data have been investigated [16][17][18][19], no methods are available for estimating synchronous responses of multiple neurons to a stimulus in a single trial when these neurons are subject to the activity that is largely unrelated to the stimulus.
In the analysis of single-trial data, it is critical to consider dependency of the current activity of neurons on the past history of their activity. A neuron undergoes an inactivation period known as a refractory period after it generates an action potential. Therefore, a model neuron significantly improves its goodness-of-fit to data if it captures this biophysical property [20,21]. In addition, estimating the dependency of the current activity level of a neuron on past spiking history of another neuron allows us to construct effective connectivity of the network within an observed set of neurons [22,23]. Including the spike-history effects in the models of synchronous ensemble activity is thus an important topic, and investigated also in [24] in the framework of a continuous-time point process theory.
In this study we construct a method for simultaneously estimating the stimulus and spikehistory effects on ensemble spiking activity when the activity of these neurons is dominated by ongoing activity. For this goal, we extend the previously developed state-space model of neuronal interactions: We model the ongoing activity, i.e., time-varying spike rates and interactions, of neurons as a latent process, and include the stimulus and spike-history effects on the activity as exogenous inputs to the latent process. We develop an expectation-maximization (EM) algorithm for this model, which efficiently combines construction of a posterior density of the latent process and estimation of the parameters for stimulus and spike-history effects. The method is tested using simulated spiking activity of 3 neurons with known underlying architecture. We provide an approximation method for determining inclusion of these exogenous inputs in the model and a surrogate method to test significance of the estimated parameters.

Methods
In this study, we analyze spike sequences simultaneously obtained from N neurons. From these spike sequences, we construct binary spike patterns at discrete time steps by dividing the sequences into disjoint time bins with an equal width of ∆ ms (in total, T bins). The width ∆ determines a permissible range of synchronous activity of neurons in this analysis. We let X t i be a binary variable of the i-th neuron (i = 1, 2, . . . , N ) in the t-th time bin (t = 1, 2, . . . , T ).
Here a time bin containing '1' indicates that one or more spikes exist in the time bin whereas '0' indicates that no spike exists in the time bin. The binary pattern of N neurons at time bin t is denoted as X t = [X t 1 , X t 2 , . . . , X t N ] ′ . The prime indicates the transposition operation to the vector. The entire observation of the discretized ensemble spiking activity is represented as

The model of time-varying simultaneous interactions of neurons
We analyze the ensemble spike patterns using time-dependent formulation of a joint probability mass function for binary random variables. Let x i be a binary variable, namely x i = {0, 1}. The joint probability mass function of N -tuple binary variables, x = [x 1 , x 2 , . . . , x N ], at time bin t (t = 1, 2, . . . , T ) can be written in an exponential form as Here θ t = [θ t 1 , θ t 2 , . . . , θ t 12 , θ t 13 , . . . , θ t 1···N ] ′ summarizes the time-dependent canonical parameters of the exponential family distribution. The canonical parameters for the interaction terms, e.g., θ t ij (i, j = 1, . . . , N ), represent time-dependent instantaneous interactions at time bin t among the neurons denoted in its subscript. ψ(θ t ) is a log normalization parameter to satisfy p(x|θ t ) = 1.
Using a feature vector that captures simultaneous spiking activities of subsets of the neurons, the probability mass function (Eq. 1) is compactly written as p(x|θ t ) = exp θ ′ t f (x) − ψ(θ t ) . The expected occurrence rates of simultaneous spikes of multiple neurons is given by a vector where expectation is performed using p(x|θ t ).
Eq. 1 specifies the probabilities of all 2 N spike patterns by using 2 N − 1 parameters. One reasonable approach to reduce the number of parameters is to select and fix interesting features in the spiking activity, and construct a probability model that maximizes entropy. For example, maximization of entropy of the spike patterns given the low-order features, yields a spin-glass model that is similar to Eq. 1, but does not include interactions higher than the second order. While it is important to explore a characteristic feature vector to neuronal ensembles, here we note that the method developed in this study does not depend on the choice of the vector, f . Below, we denote d as the number of elements in the vector, f .
Given the observed ensemble spiking activity X 1:T , the likelihood function of θ 1:T = [θ 1 , θ 2 , . . . , θ T ] is given as assuming conditional independence across the time bins. Eq.2 constitutes an observation equation of our state-space model.

Inclusion of stimulus and spike-history effects in the state model
The main focus of attention in this study is modeling of a process for the time-dependent canonical parameters, θ t , in Eq. 1. We model their evolution as a first-order auto-regressive (AR) model. The effects of the stimulus and spike history are included as exogenous inputs to the AR model (an ARX model). In its full expression, the state model is written as for t = 2, . . . , T . Here the matrix F (d × d matrix) is the first order auto-regressive parameter. ξ t (d × 1 matrix) is a random vector independently drawn from a zero-mean multivariate normal distribution with covariance matrix Q (d × d matrix) at each time bin. The state process starts with an initial value θ 1 that follows a normal distribution with mean µ (d × 1 matrix) and covariance matrix Σ (d × d matrix), namely θ 1 ∼ N (µ, Σ). Below, we describe details of the exogenous terms. The second term represents responses to external signals, or stimuli, S t , which are observed concurrently with the spike sequences. The vector S t is a column vector of n s external signals at time bin t. The each element is the value of an external signal at time bin t. If an external signal is represented as a sequence of discrete events, we denote the corresponding element of S t by '1' if an event occurred within time bin t and '0' otherwise. Multiplying S t by the matrix G (d × n s matrix) produces weighted linear combinations of the external signals at time bin t.
The third term represents the effects of spiking activity during the previous p time bins, X t−i (i = 1, . . . , p), on the current activity. The matrix H i (d × N matrix) represents the spikehistory effects of spiking activity in the previous time bin t − i on the state at time bin t. The spike-history effects are collectively denoted as Eq. 3 constitutes a prior density of the latent process in our state-space model. We denote the set of parameters in the prior distribution, called hyper-parameters, as w ≡ [F, G, H, Q, µ, Σ]. In this study, we refer to w as a parameter. In addition, we simplify Eq. 3 as θ t = Fθ t−1 + Uu t + ξ t , where u t is a single column vector constructed by stacking the stimulus vector and spike-history vectors at time bin t in a row, i.e., u t = [S t ; X t−1 ; X t−2 ; . . . ; X t−p ]. Similarly, we define a matrix U as U = [G, H]. With this simplification, the prior density defined in Eq. 3 is written as p(θ 1:T |w) = p(θ 1 |µ, Σ) T t=2 p(θ t |θ t−1 , F, U, Q), where the transition probability, p(θ t |θ t−1 , F, U, Q), is given as a normal distribution with mean Fθ t−1 + Uu t and covariance matrix Q.

Estimation of stimulus responses and spike-history effects
We estimate the parameter, w, based on the principle of maximizing a (log) marginal likelihood function. Namely, we select the parameter that maximizes For this goal, we use the expectation-maximization (EM) algorithm [25][26][27]. In this method, we iteratively obtain the optimal parameter w * that maximizes the lower bound of the above log marginal likelihood. This alternative function, known as the expected complete data log-likelihood (a.k.a., q-function), is computed as The expectation, E[ |X 1:T , w], in Eq. 5 is performed using the smoother posterior density of the state obtained by a nominal parameter, w, namely In particular, Eq. 5 can be computed using the following expected values by the posterior density: The smoother mean . These values are obtained using the approximate recursive Bayesian filtering/smoothing algorithm developed in [9,10] (See Appendix A and Eqs. A.5, A.6, and A.7 therein).
In the EM-algorithm, we obtain the parameter that maximizes the q-function by alternating the expectation (E) and maximization (M) steps. In the E-step, we obtain the above expected values in Eq. 5 by the approximate recursive Bayesian filtering/smoothing algorithm using a fixed w (Appendix A). In the M-step, we obtain the parameter, w * , that maximizes Eq. 5. The resulting w * is then used in the next E-step. Below, we derive an algorithm for optimizing the parameter at the M-step.
For the state model that includes the auto-regressive parameter and stimulus and/or spikehistory effects, these parameters are estimated simultaneously. From ∂ ∂F * q (w * |w) = 0, we obtain Here, θ t|T , W t|T , and W t−1,t|T are the smoother mean and covariance, and the lag-one covariance matrix given by Eqs. A.5, A.6, and A.7, respectively. Similarly, from ∂ ∂U * q (w * |w) = 0, we obtain Hence, the simultaneous update rule for F * and U * is given as Here the inverse matrix on the r.h.s. is obtained by using the blockwise inversion formula: The covariance matrix, Q, can be optimized separately. From ∂ ∂Q * q (w * |w) = 0, the update rule of Q is obtained as Finally, the mean of the initial distribution is updated with µ * = θ 1|T from ∂ ∂µ * q (w * |w) = 0. The covariance matrix Σ for the initial parameters is fixed in this optimization.

Simulation of a network of 3 neurons
In order to test the method, we simulate spiking activity of 3 neurons that possess specific characteristics in spike generation and connectivity as follows (See Fig. 1A). (1) The instantaneous firing rate of each neuron model depends on its own spike history in order to reproduce refractoriness in neuronal spiking activity. To achieve this, we adopt a renewal point process model whose instantaneous inter-spike interval (ISI) distribution is given by an inverse Gaussian distribution as a model of the stochastic spiking activity. (2) The firing rate of each neuron model varies across time in order to reproduce the ongoing activity. For that purpose, spike times of each neuron are generated from the renewal process by adding inhomogeneity to the underlying rate using the time-rescaling method described in [20]. The underlying rate of the inhomogeneous renewal point process model is modulated using a sinusoidal function (frequency: 1 Hz; mean and amplitude: 30 spikes/s). This rate modulation is common to the 3 neurons. (3) The neurons are activated by externally triggered stimulus inputs. To realize the stimulus responses, we deterministically induce spikes at predetermined timings of the stimuli. We consider two stimuli, one (Stimulus 1) that induces a spike in Neuron 1, and the other (Stimulus 2) that induces synchronous spikes in Neuron 2 and Neuron 3. The timings of external stimuli are not related to the sinusoidal time-varying rate, but randomly selected in the observation period (On average each stimulus happens once in 1 second). (4) There is feedforward circuitry in the network. We assume that Neuron 1 makes excitatory synaptic contacts to Neuron 2 and Neuron 3. To realize this, 5ms after a spike occurs in Neuron 1, we induce simultaneous spikes in Neuron 2 and Neuron 3 with a probability 0.5.
We simulate spike sequences with a length of 30 seconds using 1 ms resolution for numerical time steps (An example of a short period (1 s) is shown in Fig. 1B). Figure 1C displays the instantaneous spike-rates (conditional intensity functions of point processes) of Neuron 1 (Top, red line) and Neuron 2 & 3 (Bottom, green and blue lines) underlying the spiking activity in Fig. 1B. The black lines indicate sinusoidal rate modulation common to all neurons. In addition, spikes are induced in Neuron 1 at the onsets of Stimulus 1 (magenta triangles). Similarly, simultaneous spikes of Neuron 2 and Neuron 3 are generated at the onsets of Stimulus 2 (cyan triangles). In the traces of instantaneous spike-rates in Fig. 1C, instantaneous increases caused by the stimuli and synaptic inputs are not displayed. The instantaneous spike-rate of a neuron is reset to zero whenever a spike is induced in that neuron.

Selection of a state model
We analyze the simulated ensemble activity by the proposed state-space model. For this goal, we first construct binary spike patterns, X 1:T , from the simulated spike sequences of 30 seconds (Note: spike times are recorded in 1 ms resolution) by discretizing them using disjoint time bins with 2 ms width. We then apply state-space models to the binary data. The observation model used here contains interactions up to the second order (a pairwise interaction model): For the state model, we consider 5 different models that include a set of different components in Eq. 3. We select a model based on the framework of model selection in order to avoid over-fitting of a model to the data. Details of each state model are described as follows. The first state model assumes F = I, where I is an identity matrix, and does not include any of exogenous inputs. In this model, we optimize only the covariance matrix Q. The first model is denoted as [Q]. The second state model, denoted as [Q, F], is the first-order auto-regressive model. In this model, we optimize both the covariance matrix Q and the auto-regressive parameter F. The third model, denoted as [Q, F, G], additionally includes the stimulus term as exogenous inputs (Stimulus 1 and Stimulus 2). Both the matrix F and G are optimized simultaneously in addition to Q. The fourth model includes both stimulus and spike-history terms. In this model, the state model includes the history of spiking activity up to the last 6 time bins (p = 6). All parameters F, G, and H are optimized simultaneously in addition to Q. This model is denoted as [Q, F, G, H6]. The structure of the fifth model is the same as the fourth model, but contains the history of spiking activity up to the last 12 time bins (p = 12). The last model is denoted as [Q, F, G, H12].
In order to select the most predictive model among them, we select the state-space model that minimizes the Akaike (Bayesian) information criterion (AIC) [28], where w * is the optimized parameter in the Methods section. The (marginal) likelihood function in Eq. 12 is obtained by a log-quadratic approximation, i.e, the Laplace method [10] (See Appendix B for the complete equation). Figure 2B displays decreases in AICs (∆AIC) of the last four models from the AIC of the first model, [Q]. The larger the ∆AIC is, the better the state-space model is expected to predict unseen data. For these data sets, inclusion of exogenous inputs, in particular the spike history, significantly decreases the AIC. From this result, we select the state model that includes the stimulus response term and the spike-history terms up to the previous 6 time bins.

Parameter estimation
We now look at the estimated parameters of the model selected by the AIC, namely [Q, F, G, H6]. Due to the limitation in the space, we do not display the estimated dynamics of the canonical parameters, θ t , by the recursive Bayesian method (See [9,10] for the detailed analysis on dynamics of θ t by this method). The estimated parameters, G and H, are summarized in Fig. 3.
Here, in order to test the significance of the estimated parameters, we construct confidence bounds of the estimates, using a surrogate method. In this approach, we apply the same statespace model, [Q, F, G, H6], to the surrogate data set for the exogenous inputs. In the surrogate data set, the onset times of external signals (Stimulus 1 and Stimulus 2) are randomized in the    observation period. Similarly, we randomly select p = 6 bins from the past spiking activity to obtain surrogate spike history, instead of selecting the last consecutive 6 bins from time bin t. Thus the estimated parameters, G and H, are not related to the structure specified in the Section 4.1. We repeatedly applied the state-space model to the surrogate data (1000 times) to obtain the 95% confidence bound for the parameter estimation (vertical ticks in Fig. 3A, C and D). Figure 3A displays effects of the two stimuli, G, on the respective elements in θ t . First, Stimulus 1 significantly increases θ   23 by Stimulus 2 indicates that the presence of Stimulus 2 induces excess simultaneous spikes in Neuron 2 and Neuron 3 more often than the chance coincidence expected for the two neurons.
The spike-history effects are summarized as a summed matrix, p i=1 H i , shown in Fig. 3B. Two major effects are observed. First, the spike history of Neuron i significantly decreases θ (t) i (i = 1, 2, 3) (See diagonal of the first 3×3 matrix in Fig. 3B). Figure 3C displays the contribution of a spike in Neuron i during the previous p time bins to the parameter θ (t) i . These components primarily, albeit not exclusively, capture the renewal property of the simulated neuron models. Second, the spike history of Neuron 1 increases θ (t) 2 , θ (t) 3 , and θ (t) 23 (See the first column in Fig. 3B), indicating that spike interactions from Neuron 1 to Neuron 2 & 3. In particular, the increase in θ (t) 23 due to the spikes in Neuron 1 during previous 1-3 time bins (Fig. 3D) indicates that the inputs from Neuron 1 induces excess synchronous spikes in the other two neurons with approximately 2-6 ms delay.

Conclusion
We developed a parametric method for estimating stimulus responses and spike-history effects on the simultaneous spiking activity of multiple neurons when the ensemble themselves exhibit ongoing activity. The method was tested by simulated multiple neuronal spiking activity with known underlying architecture. We provided two methods to corroborate the fitted models. First, based on the result in the preceding paper, we provided an approximate equation for the log marginal likelihood (see Appendix B), which was used to select the most predictive statespace model. Second, we provided a method for obtaining confidence bounds of the estimated parameters based on a surrogate approach.
Example spike sequences simulated in this study are overly simplified. Therefore, the method needs be tested using real neuronal spike data, e.g., from cultured neurons whose underlying circuit is identified by electrophysiological studies. In practical applications, it is recommended to utilize basis functions such as raised cosine bumps used in [23] in the exogenous terms in order to capture the stimulus and spike-history effects with a fewer parameters. In addition, an appropriate bin size must be selected in order to obtain a meaningful result in the analysis of real data. Since the bin size determines a permissible range of synchronous activity, a physiological interpretation of the result depends on the choice of the bin size. It is thus recommended to present results based on multiple different bin sizes in order to confirm a specific hypothesis in a study as shown in [10,29]. Methods to overcome an artifact due to the disjoint binning are discussed in [30][31][32][33][34]. In future, inclusion of such advanced methods will allow us to detect near-synchronous responses without sacrificing temporal resolution of the analysis.
Given that applicability of the method is confirmed in real data, the proposed method is useful to investigate how ensemble activity of multiple neurons in a local circuit changes configurations of their simultaneous responses (synchronous responses) to different stimuli applied to an animal. Further, it would be interesting to see different effects of the same stimulus on the ensemble activity when an animal undergoes different cortical states.