Rhythmic activities in neurons and neuronal networks are thought to contribute to information transfer and processing in the brain. Specific frequency bands in rhythmic activities contribute to or disrupt neural information transfer from one neuron to another or from one brain region to another (Mallet et al., 2008; Sohal et al., 2009; Jackson et al., 2011; Liebe et al., 2012; McGinn and Valiante, 2014; Bastos et al., 2015). Synchronization, desynchronization, and phase-locking in rhythmic activities might play an important role in neural information representation and processing (Womelsdorf et al., 2006, 2007; van Elswijk et al., 2010; van Wingerden et al., 2010; Vinck et al., 2010; Park et al., 2013). To understand the mechanisms underpinning such functions, it is necessary to mathematically describe and analyze the properties of the underlying nonlinear dynamics (Hoppensteadt and Izhikevich, 1997). A mathematical description of this sort can be achieved with a detailed model based on biophysical mechanisms such as the kinetics of ion channels responsible for a given change in neuronal membrane potential (Hodgkin and Huxley, 1952). Such a model generally consists of a considerable number of variables and parameters. While such a detailed model enables us to directly compare its behavior with experimental data recorded under various experimental conditions, it is usually too complicated to extract essential properties of the dynamic system. As long as the detailed model possesses the nature of a nonlinear oscillator, however, its complicated dynamics can be reduced to those described by a single variable, the phase (Ermentrout and Kopell, 1984; Kuramoto, 1984; Hansel et al., 1995). Although this reduction makes impossible a quantitative comparison of the behavior of the reduced phase dynamics with experimentally observed phenomena, it provides a mathematically tractable framework through which to understand the essential properties of the dynamic system (see Figure 1a). Thus, most previous studies have first described the detailed dynamics of the model before reducing it to the simpler phase dynamics in order to analyze its dynamic properties as a rhythmic system.

However, it is quite difficult to exactly construct a detailed model based on experimental data due to nonlinearity, model selection, and parameter tuning. Therefore, it is best to obtain the reduced phase dynamics directly from experimental data if possible (Figure 1b; Tokuda et al., 2007; Kralemann et al., 2008; Cadieu and Koepsell, 2010; Stankovski et al., 2012). Ota and Aoyagi (2014) proposed a method by which one can extract the phase dynamics directly from the time-series data of coupled nonlinear oscillators using Bayesian estimation. In the proposed method, it is assumed that time-series data of all units composing a coupled oscillatory system are available. On applying the method to neural activity, however, simultaneous recoding of membrane potentials and time-series data of neuronal states from multi-neurons is still difficult, and only multi-neuronal spikes that represent a specific timing or phase in oscillatory activity can be obtained. Therefore, it remains to be established whether the proposed Bayesian approach is an effective means of extracting the phase dynamics from multi-neuronal spike data. In the present study, we addressed this issue by applying the Bayesian approach to the simulated spike data of a conductance-based neuronal network model (Figure 1c). Using the simulated spike data, we were able to compare the phase dynamics estimated directly from the multi-neuronal spike data with those theoretically derived from the detailed model.

We first conducted a numerical simulation of the computational neuronal network model in order to generate multi-neuronal spike data. We also theoretically derived the phase dynamics from the detailed model of the neuronal network with the phase reduction technique (see Phase response analysis). We then applied the Bayesian estimation to the simulated spike data and confirmed the validity of the estimated results in comparison with the reduced theoretical model (Ota and Aoyagi, 2014).

Fujita et al. (2012) reported their network to exhibit two distinct states for different parameter sets: global synchrony (case A) and two clusters corresponding to in-phase and antiphase synchrony (case B). Figures 2, 3 display the firing patterns of the network model for various conditions for cases A and B, respectively. Figure 2A is the raster display of the neurons in the network under the conditions of case A for different conditions of σ_I and σ_N. When all the neurons exhibited the same firing rate or period (σ_I = 0), the network globally synchronized irrespective of the presence of noise. In contrast, when the firing periods varied (σ_I = 0.1), the neurons no longer showed global synchronization. As shown in Figure 2B, firing periods were distributed on the basis of different values of σ_I (σ_I = 0.0: 29.03 ± 0.00 [mean ± standard deviation] ms; 0.1: 31.10 ± 2.32 ms; 0.2: 32.32 ± 5.32 ms; 0.3: 34.37 ± 9.80 ms). Figure 2C illustrates the changes in average coherence measure with or without noise as a function of bin size (σ_I = 0). Although the precision of synchronization reduced with an increase in noise intensity, the neurons synchronized even with increased noise intensity. The network activity patterns for case B are displayed in Figure 3A. When the firing periods equalized, neuronal spikes phase-locked with finite phases irrespective of noise. However, such phase-locking was lost when the firing periods varied as shown in Figure 3B (σ_I = 0.0: 29.40 ± 0.05 ms; 0.1: 30.08 ± 2.49 ms; 0.2: 31.55 ± 5.95 ms; 0.3: 34.37 ± 12.07 ms). The changes in the average coherence measure increased linearly, suggesting that the spikes were uniformly distributed (Figure 3C). As shown in Figures 2B, 3B, increased variance in the applied current changed the neuronal firing periods. Since a change in the firing period can yield different stable solutions, we confirmed the stability of the phase differences by applying the phase response analysis to different firing periods (see Materials and Methods). Figures 4A,B illustrate the odd parts of the interaction functions for cases A and B, respectively. For case A, the stable phase difference was only 0, indicating that only in-phase synchrony is stable for the examined firing periods. For case B, both 0 and π were the stable phase differences, suggesting that in-phase and antiphase synchrony are bistable. In both cases, stable phase differences were not affected by changes in firing period.

We applied the Bayesian approach to the spike data to estimate the interaction functions between the neuron pairs. Figure 5 shows the estimated interaction functions between a connected pair and an unconnected pair in case A. Figure 5A shows the estimated interaction function from neuron #32 to neuron #1 for several different numbers of spike cycles. Because there was a synaptic connection in that direction between the neurons, the interaction function theoretically obtained from the detailed model is indicated in Figure 5A by a dashed line. As more data were used, the estimated function took a shape more analogous to the theoretically derived function. In contrast, the network had no synapse to transmit a signal from neuron #17 to neuron #1. There was thus no interaction function in that direction. The estimated function thus took on a flat shape when 1,000 cycles of spike data were used. This confirms the estimation of the unconnected synapse to have been correct (Figure 5B). The shape of the interaction function determines the stability of the network states. Because the odd part of the function, Γ_odd, describes the stable states, we compared Γ_odd of the estimated function with that of the theoretically obtained function in Figure 5C. Although the scales were different, the estimated function yielded the same stable phase difference (Δϕ = 0) as the theoretical function. Thus, if the data of a sufficient number of spike cycles are available, the Bayesian approach can successfully estimate the interaction function with the parameters of case A. Estimated functions were evaluated by L²-distance between theoretically derived and experimentally derived interaction functions for various parameters (see section Materials and Methods). First, we determined the dependence of the distance on the magnitude of σ_I for σ_N = 0.4 (Figure 5D). All distances except for σ_I = 0.0 were improved as the data lengths (cycles) became larger. When σ_I = 0.0, the distance was large for all data lengths, which indicates that the proposed method did not perform well in this case. This poor estimation was independent of whether a neuron pair had direct connections (Figure S1). Figure 5E shows the dependence of the distance on the magnitude of σ_N for σ_I = 0.1. It is suggested that the distance was not significantly affected by σ_N when compared to σ_I. In this case, the distance increased as σ_N increased (red, 1,000 cycles).

Figure 6 is similar to Figure 5, except that it describes case B instead of case A. Figures 6A,B illustrate the estimated interaction functions of a connected neuronal pair (from #28 to #1) and an unconnected pair (from #12 to #1), respectively. In addition, so long as data with a sufficient number of spike cycles were used in case B, the estimated interaction functions took on the expected shape (i.e., analogous to the shape of the theoretical functions for a connected pair or a flat shape for an unconnected pair). Similar to case A, although the estimated Γ_odd was different from the theoretical one, the stabilities of the two solutions for phase differences 0 and π were unchanged (Figure 6C). Similar to Figures 5D,E,6D,E show the dependence of the L²-distances on σ_I and σ_N, respectively. As shown in Figure 6D, except for the case with 500 cycles, the distances for all values of σ_I were smaller than those in case A (Figure 5D). Even when σ_I = 0.0, the phase differences were distributed (Figures 3A,C) due to random connections, which enabled us to collect information regarding the phase relations utilized for reconstruction of the interaction functions (data not shown). The distance increased as σ_N became larger in case B as well.

The effect of partial observation on estimation performance was investigated. In the above analyses, we utilized spike data from all of the neurons in the model. Here, the results pertain to spikes recorded from a subset of modeled neurons. Figures 7A,B show the L²-distance averaged over the sampled neurons as a function of the number of sampled neurons for cases A and B, respectively. The estimation performance deteriorated with decreasing numbers of sampled neurons. Figures 7C,D illustrate how the averaged distances depended on σ_I when the numbers of sampled neurons were 32 and 16, respectively. Although the distance depends highly on both the number of cycles and σ_I, estimation performance for the case with 32 neurons was better than that of the case with 16 neurons for all values of σ_I when sufficient data were available (black, 5,000 cycles). Similarly, Figures 7E,F illustrate the dependence on σ_N in the cases with 32 and 16 neurons, respectively. In these analyses, the average distances for 32 neurons were slightly smaller than those for 16 neurons when we considered 5,000 cycles of spike data. In this case, the estimation result did not depend on whether a neuron pair had a direct connection. Subsequently, we examined whether the dependence of the estimation performance on partial sampling was the case for larger-sized networks. We increased the number of neurons in the network to 128 and 256. With an increase in the network size, the number of connections was increased whereas the synaptic strengths were decreased (see section Materials and Methods for details). Under this condition, the firing periods in the network of 128 neurons exhibited 30.94 ± 2.68 ms for case A and 29.87 ± 2.91 ms for case B. Similarly, those in the network of 256 neurons were 30.95 ± 2.27 ms for case A and 29.89 ± 2.42 ms for case B. Figures 7G,H demonstrate how the L²-distances depended on the network size for cases A and B, respectively. Similar to the case of the 64-neuron network, the estimation deteriorated as the number of sampled neurons (or the sample rate) was reduced. Furthermore, as the network size was increased, the L²-distances scaled by the synaptic strengths were increased, suggesting that the estimation performance deteriorated. This is because the synaptic strengths were decreased for the larger-size networks and it was harder for the proposed method to detect changes in phase differences for such a case.

Finally, we examined whether it is possible to infer synaptic connections based on the estimated interaction functions. We hypothesized that the summed powers of the Fourier coefficients of the estimated interaction function ∑m=1Mi{(aij(m))2+(bij(m))2} would be a good criterion for the presence of a synaptic connection. We calculated the normalized summed powers of the function P_1j for 63 pairs of neurons (from neuron #j to neuron #1, j = 2, 3,…, 64). We then applied the Otsu method to these values to classify them into two clusters based on the threshold. When P_1j was larger than the threshold, it was inferred that a synaptic connection was present from neuron #j to #1. Figure 8A shows the distribution of P_1j for case A. The top and bottom panels of Figure 8A show the results of fitting for 100- and 5,000-cycle spike data, respectively. The summed powers of neuron pairs with (red) and without (blue) actual connections were merged for 100-cycle data, but were clearly separated for 5,000-cycle data. In Figure 8B, the distributions of the summed powers for case B are shown. These data suggest that 100-cycle spike data were not sufficient for the inference of synaptic connections, while 5,000-cycle spike data were useful for this purpose. The correlation between actual and inferred connections was evaluated using MCC. In Figure 8C, MCC values are shown as functions of cycle number. Although the inferred connections did not correlate with actual connections for 100-cycle data, our method worked perfectly for 1,000-cycle data.

Constructing a detailed dynamic model based on measured data is often difficult and involves the selection of an appropriate model, adjustment for many unknown parameters, and insufficient information regarding measured data used for modeling. Furthermore, incorrect modeling can lead to misunderstanding of the properties of the dynamic system. For cases wherein the dynamic system consisted of nonlinear oscillators, Ota and Aoyagi (2014) have proposed a Bayesian method through which the reduced dynamics of a complicated dynamic system can be derived directly based on measured data (Tokuda et al., 2007; Kralemann et al., 2008; Cadieu and Koepsell, 2010; Stankovski et al., 2012). This approach is based on the fact that such oscillator systems can be represented by the reduced dynamics of phase oscillators (Ermentrout and Kopell, 1984; Kuramoto, 1984; Hansel et al., 1995). The method requires time-series data to construct the phase dynamics because the time-series data contain detailed phase information of oscillatory units. This method would be a powerful tool for the analysis of experiments recoding time-series data such as that of the electroencephalogram. However, it is still difficult to record time-series data, such as membrane potentials, from multiple neurons to investigate neural activity.

To apply the Bayesian method to spike data, it is necessary to interpolate phase values between spikes because spike timing only provides a specific phase value. Because no detailed information regarding the interspike phase values was available, we simply employed linear interpolation. This simple approximation method can be used to extract only a limited amount of information regarding the dynamics of each cycle. In addition, the neuronal network model used in the present study was rather complex and included 10 different ionic currents. Nevertheless, we were also able to use the Bayesian method to estimate the interaction function for the different conditions leading to the different network states (Figures 5, 6). Additionally, using the estimated interaction functions to infer the synaptic connections worked perfectly, although this was the case only when spike data with high cycle numbers were used (Figure 8). The method performed poorly when the network exhibited global synchronization (Figure 5 and Figure S1), when low cycle number spike data were used (Figures 5–8), when the fraction of recorded neurons was small (Figure 7), and when the synaptic strengths were weak (Figure 7). Thus, even though the data available to extract the phase dynamics consisted only of spikes, the effectiveness of the method was confirmed for several different conditions.

Estimation performance depended on firing period variability, noise intensity, and the number of sampled neurons. Although a higher degree of noise tended to degrade the estimation performance, noise intensity had a limited impact on estimation performance relative to the other factors when sufficient spike data (>5,000 cycles) were available. To examine the in vivo firing rate of GP neurons, which ranges from 30 to 50 spikes/s, a much greater input than the threshold current was applied. This resulted in a modest impact of noise on the spiking patterns and estimation performance.

Variation in firing period greatly influenced estimation performance (Figure 5D vs. Figure 5E, Figure 6D vs. Figure 6E, Figure 7C vs. Figure 7E, and Figure 7D vs. Figure 7F). Wide sampling of data on phase relations is required to improve Bayesian estimation. In the case of an identical firing period (σ_I = 0), the neurons exhibited in-phase synchronization (Figure 2A) or phase-locking with finite phase differences (Figure 3A). This limited sampling to specific phase differences. As a result, the estimation performance for the former case was the lowest among the examined conditions (Figures 5D, 7C,D). This is because when the neurons were globally synchronized, all phase differences were kept at 0. This led to biased sampling of phase relations, which is supported by the fact the confidence intervals did not converge, even for the larger cycles (Figures S1A,B). In contrast, when the variability of the firing period was not zero, but was small, it was possible to collect phase information for neuron pairs evenly over a cycle due to the slightly different firing periods. For the latter case, the phase differences for many neuron pairs were kept fixed, although those for some pairs slightly fluctuated due to random connections. This prevented us from sampling only limited phase differences. Thus, except in the special case of global synchrony, the proposed method performed well when sufficiently large and stationary spike data were available.

Partial observation is essential because it is impossible in principle to observe the activity of all neurons in the brain. Indeed, as the number of sampled neurons decreased, estimation performance deteriorated (Figures 7A,B). This is partially because the probability of sampling connected neuron pairs was decreased. If spikes of presynaptic neurons that contribute to the postsynaptic dynamics are missing (i.e., many Γ_ij(Δϕ_ij)s in Equation (2) are missing), the method would be forced to make them up, and consequently to make a wrong estimation. However, in the model investigated here, utilization of more data could serve as partial compensation. The present network model was a very small network compared to real neuronal networks. When applying the proposed method to a real neuronal network, the fraction of observable neurons should be very small. The analyses for the networks with different sizes also exhibited the same tendency that the estimation performance would deteriorate with a decrease in the fraction of observable neurons (Figures 7G,H) although we should note that the result was partly due to weakening of synaptic strengths. Thus, considering the analysis of spike data from a real neuronal network that consists of thousands or tens of thousands of neurons, the present method would be not practical and some improvements should be required. Possibly, as discussed above, the estimation performance might be improved if we can use the knowledge on anatomical connectivity to avoid randomly sample neurons to be recorded.

In the network model, neurons were connected by synapses of the same strength (synaptic conductance) or were unconnected, allowing the successful estimation of the interaction function and inference of synaptic connections so long as spike data from all neurons were used (Figure 8). Such binary connectivity potentially oversimplifies synaptic connectivity, given that synaptic strength is known to vary (Song et al., 2005). The magnitude of the interaction function is in proportion to the corresponding synaptic conductance. When using differing synaptic connection strengths, the interaction function of a weak synaptic connection could be estimated, although a large amount of spike data is required. In practice, however, it should be difficult to estimate such a weak interaction function based on limited spike data. Furthermore, it is significantly harder to make inferences when using differing synaptic connection strengths because of the inherent difficulty in separating the distributions of the summed powers. In particular, it would be quite difficult to discriminate a neuron pair with a very weak synapse from a neuron pair without synapses. Distinguishing between weak and unconnected synapses as accurately as possible would be more useful when introducing group lasso regularization to the Bayesian method.

KS conducted the data analyses and wrote the paper. TA designed the study and edited the manuscript. KK designed the study, conducted the simulations, and wrote the paper.