 Methodology
 Open Access
 Published:
Extracting singletrial neural interaction using latent dynamical systems model
Molecular Brain volume 14, Article number: 32 (2021)
Abstract
In systems neuroscience, advances in simultaneous recording technology have helped reveal the population dynamics that underlie the complex neural correlates of animal behavior and cognitive processes. To investigate these correlates, neural interactions are typically abstracted from spike trains of pairs of neurons accumulated over the course of many trials. However, the resultant averaged values do not lead to understanding of neural computation in which the responses of populations are highly variable even under identical external conditions. Accordingly, neural interactions within the population also show strong fluctuations. In the present study, we introduce an analysis method reflecting the temporal variation of neural interactions, in which crosscorrelograms on rate estimates are applied via a latent dynamical systems model. Using this method, we were able to predict timevarying neural interactions within a single trial. In addition, the pairwise connections estimated in our analysis increased along behavioral epochs among neurons categorized within similar functional groups. Thus, our analysis method revealed that neurons in the same groups communicate more as the population gets involved in the assigned task. We also showed that the characteristics of neural interaction from our model differ from the results of a typical model employing crosscorrelation coefficients. This suggests that our model can extract nonoverlapping information about network topology, unlike the typical model.
Introduction
Information communication via spike trains of neurons in populations is a core computational process that enables many brain areas to execute their roles, which include encoding of stimuli, decisionmaking, and highlevel cognition [1]. To understand these processes, therefore, the effects of spike trains across neuronal populations must be determined according to their specific network structures [2,3,4,5,6]. For this purpose, a variety of network theoretical models have been developed for the analysis of neuronal population dynamics and the description of network topology, e.g., modules, hubs, and richclubs [7,8,9]. From analyses of pairwise correlations, a recent study showed that the functional single neuron network of the macaque monkey frontoparietal area during active behavior has a highly complex topology that includes smallworldness, hubs, and richclubs consisting of oscillatory neurons synchronized in specific frequency bands [9]. To infer neural network topology, many studies have focused on quantifying correlations in the number of spikes or in the spike times of pairwise neurons. One conventionally used quantification method is the crosscorrelogram (CCG) coefficient, which measures the cooccurrence of spikes of the pairs within a given time bin [10]. To assess temporal correlations rather than the degree of synchronous firing, alternative quantification methods can be applied, such as the correlation index [11] or spiketime tiling coefficient [12], which are derived from the temporal order of pairwise spike trains. Given sufficient data, it is also possible to make use of nonlinear methods such as mutual information [13] and transfer entropy [14] in information theory.
Measuring sets of neural responses across trials is common to estimate the firing rates of spike trains because spikes are highly variable and noisy at the singletrial level. The averaged values across trials, however, obfuscate information of population dynamics at certain times in behavioral epochs or in the trials. Since estimating neural activities in single trials is challenging because of variability and deficient sampling, several studies have exploited the dimensionality reduction method adapted in the latent variable space rather than the neuronal space to measure population activities at the singletrial level [15,16,17,18,19,20,21,22,23,24,25,26,27,28,29]. Models with independent factors underlying population activities can extract neural population states in single trials [15, 16]. Recent studies have reported analysis of singletrial population activities using linear dynamical models across time [17,18,19,20] or nonlinear dynamical models allowing switches of states at given times [21, 22]. Eventdependent dynamical systems models [22] are variant models compromising linear and nonlinear methods that allow switches at the cued times between adjacent epochs. In addition, artificial recurrent neural networks infer singletrial neural population dynamics based on the assumption that the networks can generate neural data using a machinelearning method [23]. In nonhuman primates, population activities in diverse brain regions are low dimensional in many instances [19, 30,31,32,33,34,35,36,37]; therefore, it can be assumed that population activity in single trials can be estimated via a latent dynamical systems model with much lower dimensions than the population size.
Estimations of population activity are useful for elucidating the slowvarying neural dynamics governed by shared inputs in a single trial; however, the recovery of transient changes evoked by other neurons’ spikes is limited because the temporal precision of the latent model is generally insufficient for capturing these interactions. The intrinsic timescale of spike effects is about tens of ms and below 200 ms at most [38], whereas neural populations sustain their activity in the sameshared activity space within an epoch that generally lasts longer than hundreds of ms [22]. To address this issue, it is possible to calculate CCGs of pairs of neurons by two different ways that CCGs are calculated using spike trains or rate estimates. In typical fashion, CCGs are calculated using spike trains, but they include the effects caused by shared motives such as common inputs or synchronization across the populations as well as the direct effects of the other neurons in the population. To quantify these shared terms, the CCGs can be calculated using estimates of population activity inferred by a latent dynamical systems model. The transient effects between pairs of neurons can be surmised by comparing two CCGs according to spike trains or estimates of populate activity. Estimates of population activity can then be revised in light of transient effects, and the degree to which both transient and shared effects contribute to the revised estimate can be evaluated.
In the present study, we applied the abovedescribed method to a dataset from an open database in which data were simultaneously recorded from the anterior lateral motor (ALM) cortex of mice while they were executing a twoalternativeforced choice task [29, 39,40,41,42,43]. We found that neural interactions varied significantly over time in a trial and that neurons could be classified based on the degree of other neurons’ effects in each behavioral epoch. We also showed that network structures were dependent on the taskrelated epochs, implying that interactions between neurons were not stationary but instead adjusted to the relative importance of the epochs. Compared with a model containing correlation coefficients, which are widely used to evaluate functional connectivity, our model could improve judgment of the underlying organization of transient neural activities.
Materials and methods
Application to mice ALM cortex neurons
We applied the abovedescribed analysis method to a dataset from alternative choice tasks conducted in previous studies [29, 39,40,41,42,43]. In this research, trained mice reported pole position (posterior or anterior) by licking one of two targets (left in an ipsi trial and right in a contra trial) following a delay epoch after pole presentation. After a 1.3s delay, mice responded to an auditory signal indicating the onset of the response epoch (Fig. 1a, upper panel). Extracellular spikes were recorded on the lefthemisphere ALM cortex using 32channel NeuroNexus silicon probes (19 mice) or 64channel Janelia silicon probes (20 mice). The extracellular recording traces were bandpassfiltered (300–6000 Hz) and then sorted using JRclust if events exceeded an amplitude threshold. Based on spike width, cells were classified as pyramidal cells (width > 0.5 ms) or fastspiking interneurons (width < 0.35 ms). For a more detailed explanation of the electrophysiological recording method, consult the studies cited earlier [40, 42, 43]. To assess the effects of ensemble inputs from other neurons, we analyzed only the sessions that included > 15 neurons (9 of 40 sessions).
Estimating the dynamics of ALM cortex neural activities
Neural dynamics are highly complicated and stochastic, resulting in diverse responses even in comparable cases. To resolve these difficulties, the extraction of shared neural activity (SNA) in a lowdimensional latent space has been developed to determine the dynamics of population activity as an alternative to directly assessing raw neural activities. For this purpose, we used a previously introduced model [22], i.e., an eventdependent linear dynamical systems (EDLDS) model that allows switches of behavioral epochs as follows:
The latent variable \(x\left(t\right)\) is an \({N}^{L}\)dimensional vector (\({N}^{L}\): the dimension of the latent space) updated by matrix \(W\left(s\right)\) in the sepoch (s: sample, delay, response). The term \({W}_{0}(t)\) represents random fluctuation independent of the term updated by matrix \(W\left(s\right)\) in the latent space. The rate estimate \(r(t)\) is an \({N}^{P}\)dimensional vector (\({N}^{P}\): the number of neurons in the population) calculated by the projection matrix \(V\left(s\right)\) in the sepoch. The term \({V}_{0}(t)\) denotes the residual that is not explained by the latent variables.
Temporal interaction estimated with CCGs. Since the EDLDS model captures slowvarying dynamics of the neural response, neural components that reflect temporal fluctuations must be introduced into the model. One such possible component is the interaction between neurons over a finite time scale. We assumed that a neuron was affected by all the other neurons in the population and that the rate estimate calculated in Eq. 1 could be revised by the integration of these interactions. Accordingly, we calculated CCGs between all pairs of neurons from their spike trains \({C}_{ij}(k\tau )\) as follows:
where \(j\) (\(i\)) is the index of a reference (a target) neuron affecting (affected by) the other neuron. \(M\) and \({N}_{j}^{m}\) are the number of trials and the number of spikes of the reference neuron during trial m, respectively, and \({n}_{i}^{m}\left({t}_{jl}^{m},k\right)\) denotes the number of spikes of the target neuron included in the kth time lag range \({t}_{jl}^{m} + \left(k  1\right)\tau \le t < {t}_{jl}^{m} + k\tau\) after the lth spike time \({t}_{jl}^{m}\) of the reference neuron. We found similar results with various temporal precisions from \(\tau =15 ms\) to \(\tau =45 ms\) for the transient effects; thus, we heuristically selected the temporal precision \(\tau =25 ms\) because the results of analysis using this value showed greater distinction in neural interactions according to the epoch. In a similar manner, the estimated CCG (ECCG) from the SNA of a latent model \({C}_{ij}^{E}(k\tau )\) can be represented as:
where \(\langle {r}_{i}^{m}\left({t}_{jl}^{m},k\right)\rangle\) denotes the mean rate estimate of the target neuron in the kth time lag range \({t}_{jl}^{m }+ \left(k  1\right)\tau \le t < {t}_{jl}^{m} + k\tau\) after the lth spike time \({t}_{jl}^{m}\) of the reference neuron. In line with \({n}_{i}^{m}\left({t}_{jl}^{m},k\right)/\tau\) in Eq. 3, which denotes the rate in the range \({t}_{jl}^{m} + \left(k  1\right)\tau \le t <{t}_{jl}^{m} + k\tau\) derived from the spike trains of the ith neuron, the counterpart \(\langle {r}_{i}^{m}\left({t}_{jl}^{m},k\right)\rangle\) in the ECCG (Eq. 4) is the rate in the same range derived from the rate estimate of the ith neuron (Eq. 1). The difference between the CCG and ECCG (i.e., the DCCG) \(\Delta {C}_{ij}\left(k\tau \right)\) can then be represented as follows:
Optimization of neural activities
SNA effectively describes the slowvarying population activity governed by the shared inputs. However, recovering the transient changes evoked by the spikes of the connected neurons is difficult due to the short time variability of the interaction effect compared with the temporal precision of the latent model. For increasingly accurate estimates that include these transient effects over a relatively short time scale, the rate estimate in the given time bin of the latent model can be corrected using the CCG and ECCG. Since the interaction between neurons might change across times or trials, we allowed the CCG and ECCG to depend on both variables. We calculated the CCG and ECCG from the spike trains in the neighboring 20 bins of 20 trials (we obtained similar results using different values, e.g., 15 or 30 bins or trials; data not shown). For simplicity, we used only onetime lag before and after the spike of the reference neuron (\({C}_{ij}\left(\tau \right), {C}_{ij}\left(\tau \right)\)) and then the corrected rate estimate \({r}_{i,C}^{T}(t)\) as follows:
where \({\delta }^{+}\left(\left\{{t}_{j}\right\},t\right) = 1\) or 0 (\({\delta }^{}\left(\left\{{t}_{j}\right\},t\right) = 1\) or 0) if neuron j fires (or does not fire) in the time lag range \(\left[\mathrm{t},t + \tau \right]\) (\(\left[t  \tau ,t\right]\)), T denotes the time bin index of the latent model, and \(0 \le t <{ S}_{bin}\) (\({S}_{bin}\): time bin size of the latent model). We used \({S}_{bin}\)= 67.4 ms following a previous study [22]. The rate estimate \({r}_{i}^{T}\) represents the SNA of the ith neuron in the given time bin \(T\). Because of the effects of the other neurons’ spikes, the mean of the corrected rate estimate \({r}_{i,C}^{T}(t)\) in the given time bin T becomes different from the rate estimate \({r}_{i}^{T}\). To resolve this disagreement, we normalized the corrected rate estimate to match with the SNA; the transient SNA (TSNA) \({r}_{i,TSNA}^{T}\left(t\right)\) was then represented as follows:
where \(\langle {r}_{i,C}^{T}(t)\rangle\) is the mean of \({r}_{i,C}^{T}(t)\) in the given time bin \(T\). Since the contributions of transient effects to the neural response vary depending on the time bin, we found the optimal degree of SNA and TSNA contributions by maximizing the likelihood \(L\left(\left\{{t}_{i}\right\}\right)\) of the target neuron’s spike times in a given time bin as follows:
where \(\left\{{t}_{i}\right\}\) denotes the spike trains of neuron i and \({\delta }_{t,\left\{{t}_{i}\right\}} = 1\) or 0 if \(t\in \left\{{t}_{i}\right\}\) or does not, respectively. The probability \(P(t)\) of firing at time t is calculated using the softmax selection rule \(P(t) = \left(\alpha {r}_{i,TSNA}^{T}\left(t\right) + (1  \alpha ) {r}_{i}^{T}(t)\right)\delta t\), where \(\alpha = 1 / (1 +{ e}^{Q})\). We set a discrete time step of \(\delta t={S}_{bin}/60\)(\({S}_{bin}\)= 67.4 ms: time bin size of the latent model) for all calculations in this work. The parameter \(\alpha\) denotes the degree to which TSNA contributes and it is determined by Q. The optimized neural activity (ONA)\({r}_{i,ONA}^{T}\left(t\right)\) is then represented as follows:
where \({\alpha }_{*}\) denotes the value maximizing \(L\left(\left\{{t}_{i}\right\}\right)\), i.e., \({\alpha }_{*} = 1 / (1 + {e}^{{Q}_{*}})\) and \({Q}_{*} = arg{max}_{Q} L\left(\left\{{t}_{i}\right\}\right)\).
Neural interactions estimated by Kullback–Leibler divergence
To assess how SNA varies by the transient effects of other neurons, we compared the probability distributions of the spike trains, i.e., Kullback–Leibler divergences, between the ONA and SNA (i.e., the DONA) \({D}_{T}^{KL}({P}_{T}^{ONA}\left(t\right){P}_{T}^{SNA})\) as follows:
where \({P}_{T}^{SNA}\) (\({P}_{T}^{ONA}\left(t\right)\)) denotes the probability of firing for SNA(ONA) at time t in the given time bin \(T\). To identify interactions induced by a single reference neuron’s spike train but not by all other neurons (as with DONA), we defined the pairwise Kullback–Leibler divergence between ONA and SNA (i.e., the PDONA) \({D}_{T}^{KL}\left({P}^{ONA}\left({s}_{n}\right){P}^{SNA}\left({s}_{n}\right)\right)\) as follows:
where \({s}_{n}\) is an ndimensional vector with 0 or 1 denoting spike existence or nonexistence of the target neuron i in the time step \(\delta t\) after nconsecutive spike times of the reference neuron j (\(S(n)\): the entire space including all the possible cases of \({s}_{n}\)) closest to the time bin T, whereas \({t}_{j}^{k}\) denotes the kth spike time of the reference neuron j, and \({\delta }_{1,{s}_{n}\left(k\right)}=\) 1 or 0 if the kth element of \({s}_{n}\) (\({s}_{n}(k)\)) is 1 or 0, respectively. \(P\left({s}_{n}\right)\) represents the probability of the series of the target neuron firing or not firing after the firing of the reference neuron. Hence, PDONA measures how each of the reference neuron’s nconsecutive spikes modify the probability of the target neuron firing after the moment at which the corresponding reference neuron spikes. We tested various cases of n spikes of reference neurons (from n = 5 to n = 8); we chose to use n = 6 and excluded cases in which the number of spikes by the reference neuron was < 10 in the trial (we obtained similar results for n = 5, 7, or 8; data not shown).
Logistic regression analysis for decoding lick responses
We evaluated how the neural interaction improved performance in decoding lick responses through logistic regression analysis, including spike counts and DONA (full model), as follows:
where the logodds for the probability of the lick of the left target in the ipsi trial (\({P}_{ipsi}\)) were fitted using the log of the spike counts (\(\mathrm{log}({N}_{sp})\)) and the log of DONA (\(\mathit{log}\left(DONA\right)\)). By comparing Akaike’s information criterion (AIC) between the full model and reduced model including only the spike counts, we determined the optimal model for decoding lick responses.
Crosscorrelation coefficients
To compare DONA with a typical measure of neural interaction, we calculated the crosscorrelation coefficients of neuron pairs in the given time bin. For direct comparison with DONA, we defined a normalized correlation coefficient (NCC) \({\rho }_{i}\), which assumes the mean correlation between a single neuron i and all the other neurons as follows:
where \({N}^{P}\) is the number of neurons in the population and \({\rho }_{ij}\) is the crosscorrelation coefficient between neuron i and neuron j.
Results
Population activities and CCGs
To demonstrate the neural interactions in an ALM cortex, we analyzed neural data collected from mice as they executed a delayed response task (Fig. 1a, upper panel). ALM neurons showed complex and variable activities, as shown in the example raster plot and firing rates from a single trial (Fig. 1a, middle and lower panels). For exemplification of different CCGs from the spike trains and rate estimates, we provide CCGs for three different pairs of neurons. N15 is the target neuron affected by the other neurons, i.e., the reference neurons; we show the effects of only three reference neurons, N2, N7, and N12, on the target neuron for the purpose of improving visibility. Examples of CCGs, ECCGs, and DCCGs for three reference–target pairs are shown in Fig. 1b. Both the CCGs and ECCGs were sustained over hundreds of ms as they moved farther away from the reference neuron’s spike time (time = 0). DCCGs, however, approached zero after a few bins (about tens of ms), implying that the DCCG captures the transient effect near time zero by offsetting the shared effects of both the CCG and ECCG remaining in the longer time range.
Localized neural interaction terms at specific behavioral epochs
In our analysis, DONA is a quantity with which we can judge how the probability distribution of the ONA differs from that of the SNA, which is acquired from the latent model (Eqs. 1–2). Figure 2a provides an example of the optimization of the activities of the target neuron N15 (Fig. 1b), including all the effects of the other (reference) neurons, with parameter \(\alpha\) deciding the degree of SNA and TSNA contributions (Fig. 2b). Since ONA is the optimal neural activity between the transient component from other neurons and shared components due to behavior, larger values of DONA can be interpreted as more interactions occurring with the other neurons. Neuron N15 in Fig. 2c, for example, shows more interactions with other neurons during the sensory (indicated by S) and response (indicated by R) epochs compared with the interaction observed during the delay epoch (indicated by D). For pairwise interactions, we showed PDONAs of three pairs, which represent the effects of three reference neurons N2, N7, and N12 on the target neuron N15 shown in Fig. 1b (Fig. 2d). It should be noted that the DONA of the target neuron N15 in Fig. 2c was the interaction integrated with all other neurons, although we only showed the results of the three reference neurons to improve visibility. To characterize the neural interactions of all neurons, we classified them in conformity with the epochs in which their DONAs were maximized. First, to quantify the effect over time in a trial, we averaged DONA across trials. Neurons were considered to show more interactions with other neurons in an epoch if the magnitudes of DONA in the epoch (20 bins per epoch) were larger than the DONA magnitudes in the other two epochs (Wilcoxon rank sum test, p < 0.05). In an example session shown in Fig. 3a (session no. 37), DONAs were significantly larger; therefore, they were well localized in specific epochs. Among 22 neurons, 6 (27%), 5 (23%), and 5 (23%) showed more interactions with other neurons in the sample, delay and response behavioral epochs, respectively. For the neurons of all sessions (157 neurons in 9 sessions), most (111 of 157; 70.70%) showed maximum values of DONA consistently in single epochs (Fig. 3b). This result implies that the effects of neural interactions on single neurons are rather concentrated on a moment of information processing than on being maintained throughout the entire trial. The remaining neurons showed larger values of DONA during the two epochs (26 of 157; 16.54%) or no larger values (20 of 157; 12.74%).
Tendency of increases of neural interactions toward the end of a trial
To determine the overall variation of DONA during the progress of epochs, we simply averaged DONA values for all the neurons regardless of their statistical significance in epochs (Fig. 3c). The mean DONA value increased and reached a peak soon after the response onset. Significantly more neurons showed larger DONA values in the latter epochs, in agreement with the results of mean DONA values (p < 0.05, χ^{2} test; Fig. 3d). The number of neurons with larger DONA values was 33, 50, and 68 of 157 total neurons during the sample, delay, and response epochs, respectively; thus, there were significant differences between the behavioral epochs (p < 0.05 between S and D and between D and R, both χ^{2} tests).
In the dataset used for this analysis, Hidehiko et al. [43] categorized neurons in the ALM cortex into four functional groups, namely contralateral rampingup, ipsilateral rampingup, rampingdown, and the rest based on the selectivity of the lick responses and the neurons’ ramping activities. We investigated whether the neuronal interaction term DONA was related to these functional neuron groups. According to the categorization of Hidehiko et al. [43], we placed 52, 17, 19, and 69 neurons into the contralateral rampingup, ipsilateral rampingup, rampingdown, and remaining groups, respectively. Secondly, we calculated PDONAs between neurons belonging to the same functional groups and the different functional groups separately. As was the case for DONA, pairs of neurons were considered to show larger interactions in an epoch if the magnitudes of PDONA in the epoch (20 bins per epoch) were larger than the magnitudes in the other two epochs (Wilcoxon rank sum test, p < 0.05). The proportion of pairs with larger values of PDONA in each epoch is shown Fig. 4a (filled and empty bars show pairs of neurons in the same functional group and different groups, respectively). Significantly more pairs in the same functional group showed larger PDONA values in the latter epochs in accordance with the DONA results shown in Fig. 3d (Fig. 4a, black bars). In the same functional group, the pairs of neurons showing larger PDONA values were 36, 73, and 121 of 184 total pairs during the sample, delay, and response epochs, respectively (p < 0.01 between S and D and between D and R, both χ^{2} tests). However, for pairs of neurons belonging to different functional groups, there was no significant difference in the number of pairs with larger PDONA values between epochs (Fig. 4a, empty bars): 374, 403, and 437 of 1,125 total pairs during the sample, delay, and response epochs, respectively (p = 0.199 between S and D; p = 0.138 between D and R; both χ^{2} tests). These results imply that the increase in neural interactions in the latter epoch was due to increased functional connections between neurons in the same functional group.
The number of pairs of larger PDONA values was significantly different in the sample and response epochs but not in the delay epoch (delay: p = 0.314, χ^{2} test). In the sample epoch, the proportion of pairs with larger PDONA values was significantly higher when neurons in the pairs belonged to different functional groups (p < 0.01, χ^{2} test). However, the opposite was true in the response epoch in which the proportion of pairs of neurons in the same functional group was higher (p < 0.01, χ^{2} test). We also reanalyzed the data for each group separately (Fig. 4b). The numbers of pairs of neurons in the contralateral rampingup, ipsilateral rampingup, and rampingdown group were 150, 18, and 16 of 184 total pairs of neurons. It should be noted that we did not analyze the fourth group (the rest) separately because neurons in the fourth group were just the remainders that could not be classified into any of the three groups discussed above. Although the results for the contralateral rampingup group were identical to those of all the groups shown in Fig. 4a (Fig. 4b, left panel), the PDONAs of the other groups had different distributions when compared with the overall tendency to increase along the epoch (Fig. 4b, middle and right panels). In the contralateral rampingup group, the pairs of neurons showing larger PDONA values were 26, 54, 117 of 150 total pairs during the sample, delay, and response epochs, respectively (p < 0.01 between S and D and between D and R, both χ^{2} tests). In the ipsilateral rampingup group, the pairs of neurons showing larger PDONA values were 2, 13, 4 of 18 total pairs during the sample, delay, and response epochs, respectively (p < 0.01 between S and D and between D and R, both χ^{2} tests). In the rampingdown group, the pairs of neurons showing larger PDONA values were 8, 6, 0 of 16 total pairs during the sample, delay, and response epochs, respectively (p = 0.48 between S and D; p < 0.01 between D and R; both χ^{2} tests). This suggests that neurons classified in different groups show different response patterns. In addition, the interactions between these neurons also show different tendencies according to epoch.
Enhanced decoding capability for behavioral responses due to neural interactions
We performed logistic regression analyses with the spike count \({N}_{sp}\) and DONA to test whether neural interactions enhanced the performance when decoding lick responses (Fig. 5a; see Material and methods). The probability of the response was significantly dependent (p < 0.05 in the logistic regression) on DONA in 10.53% of neurons (p < 0.05, binomial test) and on the spike count in 71.58% of neurons (p < 0.01, binomial test), as shown in Fig. 5a. We compared the full model including both DONA and spike counts, which we expected to be optimal for single neurons, with the reduced model including only spike counts. Comparing the AICs of the two models, we found that the full model was optimal in a considerable number of neurons (21.13%; Fig. 5b). Thus, it seems that DONA can convey more information than spike count alone. In Fig. 5c, d, we compare the dynamics of DONA and spike counts according to lick responses by showing a representative example of a single neuron for which the probability of the response depended significantly on both DONA and the spike counts of the neuron in the delay epoch.
Disparate trends of crosscorrelation compared with DONA
To compare the localization trends of DONA (shown in Fig. 3a–c) with an existing and typical measure of pairwise interactions, we calculated NCCs for each neuron (Fig. 6; see Materials and methods). To enable the direct comparison of these measures, the results of NCC were arranged for the neurons sorted in the same manner shown in Fig. 3a–c. If NCC values could extract similar features of network topology as DONA, we would see a resemblance between NCC and DONA values in terms of localization to epochs. However, NCCs did not show such a localization tendency and the mean NCC showed no patterns along the epochs, unlike the mean DONA (Fig. 6b, open circle). These results suggest that DONA might extract nonoverlapping features related to the neural interactions in the population that cannot be extracted by the typically used correlation measure.
Discussion
In the present study, we present a method by which to extract singletrial neural interactions using rate estimates obtained by a latent dynamical model. We applied this method to a dataset containing data from mice cortical neurons [22, 39,40,41,42,43]. The model applied in this method was able to evaluate the proportions of shared neural activities and transient effects evoked by other neurons; thus, it showed that inferred measures could trace the variation of neural interactions along epochs. Most neurons were maximally affected by other neurons in a specific epoch and, on average, neural interactions increased as an epoch progressed. To demonstrate the novelty of our method, we compared its performance to that achieved using crosscorrelation coefficients; we showed that our model could reveal hidden features of the neural network that were not accounted for by a typical method in which pairwise interactions are estimated.
Obstacles must be overcome to ensure that our method can be generally applied to neural population data in different regions of the brain when performing diverse tasks. First, a number of existing models [15,16,17,18,19,20,21,22,23,24,25,26,27,28,29] have been established by which to analyze singletrial neural population dynamics and provide different rate estimates to certain extents. It might also be more challenging to estimate neural activities in single trials for tasks that require higher computational capabilities [21], although many studies have shown that population activities could be reduced to lowdimensional shared spaces [25, 30,31,32,33,34,35,36,37]. In such tasks, latent dynamical models provide suboptimal rate estimates and the consequential inference of neural interactions will mislead the coordination of the neural network. Hence, it will be necessary to analyze how discrepancies in rate estimates change the properties of DONAs or PDONAs and how the suitability of DONAs or PDONAs can be judged when using diverse datasets with varying levels of difficulty.
Second, timescales describing SNA and TSNA will vary according to types of neuron or task [22, 38]. Furthermore, correlations between neurons are identified in accordance with the timescale of their pairwise effects [12]. Consequently, it is vital that the temporal precision of the rate estimates and the transient effects are chosen appropriately to allow confident interpretation of neural interactions. In the present study, we chose the temporal precisions heuristically so that the analysis model showed results that were more distinctive. Therefore, we cannot exclude the possibility that models with different temporal precisions might provide new information about neural interactions. Additionally, for realtime processing, such as in brain–machine interfaces or closedloop brain stimulations, temporal precisions must be automatically regulated depending on the given neural data. Thus, for the general application of our method, it will be essential to identify techniques and criteria by which to define temporal precisions according to the characteristics of the neural data.
Third, we only compared our methods to crosscorrelation coefficients. A number of other measures exist for quantifying interactions that have various uses related to the characteristics of neural data or the purposes of analysis [10,11,12,13,14]. To validate the effectiveness of DONA and PDONA as measures for inferring neural interaction, comparisons to other measures, including the degree of synchronous firing or the temporal correlation of neural pairs, will be required to identify the similarities and differences. Indeed, we should be open to the possibility that combinations of singletrial rate estimates and measures of correlations might provide more information about network structure.
We anticipate that this study will be extensible for further research into neural population dynamics. Our developed method can provide novel insights into types of neurons or neural pairs based on neural interactions, just as neurons are classified into certain types based on their functional roles or responses in specific epochs or sessions [43]. Here we applied our method only to spike data with binary values. For its application to continuous datasets, such as those containing EEG or fMRI data, further research will be required. In particular, we must be able to define the shared activities and transient effects of continuous datasets in line with the measures applied in this work.
Availability of data and materials
The datasets supporting the conclusions of this article is available in the reference [42]. Interim data processed for supporting the conclusions of this article will be made available from the authors on reasonable request.
Abbreviations
 ALM:

Anterior lateral motor
 SNA:

Shared neural activity
 EDLDS:

Eventdependent linear dynamical systems
 CCG:

Crosscorrelogram
 ECCG:

Estimated CCG
 DCCG:

Difference between CCG and ECCG
 TSNA:

Transient SNA
 ONA:

Optimized neural activity
 DONA:

Kullback–Leibler divergence between ONA and SNA
 PDONA:

Pairwise DONA
 NCC:

Normalized crosscorrelation
References
 1.
Rieke F, Warland D, Ruyter van Stevninck R, Bialek W. Spikes: Exploring the neural code. Cambridge: MIT Press; 1999.
 2.
Pouget A, Beck JM, Ma WJ, Latham PE. Probabilistic brains: knowns and unknowns. Nat Neurosci. 2013;16(9):1170–8.
 3.
Wohrer A, Humphries M, Machens C. Populationwide distributions of neural activity during perceptual decisionmaking. Prog Neurobiol. 2013;103:156–93.
 4.
Yuste R. From the neuron doctrine to neural networks. Nat Rev Neurosci. 2015;16:487–97.
 5.
Fornito A, Zalesky A, Bullmore E. Fundamentals of brain network analysis. CA: Academic Press; 2016.
 6.
Humphries MD. Dynamical networks: finding, measuring, and tracking neural population activity using network science. Netw Neurosci. 2017;1(4):324–38.
 7.
Watts DJ, Strogatz SH. Collective dynamics of ‘smallworld’ networks. Nature. 1998;393:440–2.
 8.
Bullmore E, Sporns O. Complex brain networks: graph theoretical analysis of structural and function systems. Nat Rev Neurosci. 2009;10:186–98.
 9.
Dann B, Michaels JA, Schaffelhofer S, Scherberger H. Uniting functional network topology and oscillations in the frontoparietal single unit network of behaving primates. eLife. 2016;5:e15719.
 10.
Pasquale V, Massobrio P, Bologna LL, Chiappalone M, Martinoia S. Selforganization and neuronal avalanches in networks of dissociated cortical neurons. Neuroscience. 2008;153:1354–69.
 11.
Wong RO, Meister M, Shatz CJ. Transient period of correlated bursting activity during development of the mammalian retina. Neuron. 1993;11:923–38.
 12.
Cutts CS, Eglen SJ. Detecting pairwise correlations in spike trains: an objective comparison of methods and applications to the study of retinal waves. J Neurosci. 2014;34(43):14288–303.
 13.
Singh A, Lesica NA. Incremental mutual information: A new method for characterizing the strength and dynamics of connections in neuronal circuits. PLoS Comput Biol. 2010;6:e1001035.
 14.
Thivierge JP. Scalefree and economical features of functional connectivity in neuronal networks. Phys Rev E. 2014;90:022721.
 15.
Yu BM, Cunningham JP, Santhanam G, Ryu SI, Shenoy KV, Sahani M. Gaussianprocess factor analysis for lowdimensional singletrial analysis of neural population activity. J Neurophysiol. 2009;102:614–35.
 16.
Zhao Y, Park IM. Variational latent Gaussian process for recovering singletrial dynamics from population spike trains. Neural Comput. 2017;29:1293–316.
 17.
Paninski L, Ahmadian Y, Ferreira DG, Koyama S, Rad KR, Vidne M, et al. A new look at statespace models for neural data. J Comput Neurosci. 2010;29:107–26.
 18.
Macke JH, Buesing L, Cunningham JP, Yu BM, Shenoy KV, Sahani M. Empirical models of spiking in neural populations. Adv Neural Inf Process Syst. 2011;24:1350–8.
 19.
Kao JC, Nuyujukian P, Ryu SI, Churchland MM, Cunningham JP, Shenoy KV. Singletrial dynamics of motor cortex and their applications to brainmachine interfaces. Nat Commun. 2015;6:7759.
 20.
Elsayed GF, Lara AH, Kaufman MT, Churchland MM, Cunningham JP. Reorganization between preparatory and movement population responses in motor cortex. Nat Commun. 2016;7:13239.
 21.
Morcos AS, Harvey CD. Historydependent variability in population dynamics during evidence accumulation in cortex. Nat Neurosci. 2016;19:1672–81.
 22.
Wei Z, Inagaki H, Li N, Svoboda K, Druckmann S. An orderly singletrial organization of population dynamics in premotor cortex predicts behavioral variability. Nat Commun. 2019;10:216.
 23.
Pandarinath C, O’Shea DJ, Collins J, Jozefowicz R, Stavisky SD, Kao JC, et al. Inferring singletrial neural population dynamics using sequential autoencoders. Nat Methods. 2018;15:805–15.
 24.
Ahrens MB, Li JM, Orger MB, Robson DN, Schier AF, Engert F, Portugues R. Brainwide neuronal dynamics during motor adaptation in zebrafish. Nature. 2012;485:471–7.
 25.
Kaufman MT, Churchland MM, Ryu SI, Shenoy KV. Cortical activity in the null space: permitting preparation without movement. Nat Neurosci. 2014;17:440–8.
 26.
Cunningham JP, Yu BM. Dimensionality reduction for largescale neural recordings. Nat Neurosci. 2014;17:1500–9.
 27.
Raposo D, Kaufman MT, Churchland AK. A categoryfree neural population supports evolving demands during decisionmaking. Nat Neurosci. 2014;17:1784–92.
 28.
Kiani R, Cueva CJ, Reppas JB, Newsome WT. Dynamics of neural population responses in prefrontal cortex indicate changes of mind on single trials. Curr Biol. 2014;24:1542–7.
 29.
Kaufman MT, Churchland MM, Ryu SI, Shenoy KV. Vacillation, indecision and hesitation in momentbymoment decoding of monkey motor cortex. eLife. 2015;4:e04677.
 30.
Afshar A, Santhanam G, Yu BM, Ryu SI, Sahani M, Shenoy KV. Singletrial neural correlates of arm movement preparation. Neuron. 2011;71:555–64.
 31.
Harvey CD, Coen P, Tank DW. Choicespecific sequences in parietal cortex during a virtualnavigation decision task. Nature. 2012;484:62–8.
 32.
Churchland MM, Cunningham JP, Kaufman MT, Foster JD, Nuyujukian P, Ryu SI, et al. Neural population dynamics during reaching. Nature. 2012;487:51–6.
 33.
Mante V, Sussillo D, Shenoy KV, Newsome WT. Contextdependent computation by recurrent dynamics in prefrontal cortex. Nature. 2013;503:78–84.
 34.
Sadtler P, Quick KM, Golub MD, Chase SM, Ryu SI, TylerKabara EC, et al. Neural constraints on learning. Nature. 2014;512:423–6.
 35.
Carnevale F, de Lafuente V, Romo R, Barak O, Parga N. Dynamic control of response criterion in premotor cortex during perceptual detection under temporal uncertainty. Neuron. 2015;86:1067–77.
 36.
Pandarinath C, Gilja V, Blabe CH, Nuyujukian P, Sarma AA, Sorice BL, et al. Neural population dynamics in human motor cortex during movements in people with ALS. eLife. 2015;4:e07436.
 37.
Kobak D, Brendel W, Constantinidis C, Feierstein CE, Kepecs A, Mainen ZF, et al. Demixed principal component analysis of neural population data. eLife. 2016;5:e10989.
 38.
Cirillo R, Fascianelli V, Ferrucci L, Genovesio A. Neural intrinsic timescales in the macaque dorsal premotor cortex predict the strength of spatial response coding. iScience. 2018;10:203–10.
 39.
Guo ZV, Li N, Huber D, Ophir E, Gutnisky D, Ting JT, et al. Flow of cortical activity underlying a tactile decision in mice. Neuron. 2014;81:179–94.
 40.
Li N, Chen TW, Guo ZV, Gerfen CR, Svoboda K. A motor cortex circuit for motor planning and movement. Nature. 2015;519:1–16.
 41.
Li N, Daie K, Svoboda K, Druckmann S. Robust neuronal dynamics in premotor cortex during motor planning. Nautre. 2016;532:1–25.
 42.
Guo ZV, Inagaki HK, Daie K, Druckmann S, Gerfen CR, Svoboda K. Maintenance of persistent activity in a frontal thalamocortical loop. Nature. 2017;545:181–6.
 43.
Inagaki HK, Inagaki M, Romani S, Svoboda K. Lowdimensional and monotonic preparatory activity in mouse anterior lateral motor cortex. J Neurosci. 2018;38:4163–85.
Acknowledgements
We thank members of center for neuroscience imaging research at Sungkyunkwan University for discussion and proofreading.
Funding
This research was supported by grants of the Brain Research Program through the National Research Foundation of Korea funded by the Ministry of Science and ICT (NRF2016M3C7A1904986) and by IBSR015D1.
Author information
Affiliations
Contributions
NH and JS designed research; NH, SK and JS analyzed data. All authors wrote the paper. All authors read and approved the final manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Huh, N., Kim, SP., Lee, J. et al. Extracting singletrial neural interaction using latent dynamical systems model. Mol Brain 14, 32 (2021). https://doi.org/10.1186/s13041021007407
Received:
Accepted:
Published:
Keywords
 Neural interaction
 Latent dynamical systems model
 Crosscorrelogram
 Optimized neural activity