Introduction

Understanding how information about the sensory stimulus is encoded in sequences of action potentials (spikes) of sensory neuronal population is a major focus of systems neuroscience. The recent studies focusing on exploring the functional connectivity property of the brain network (Gupta et al. 2011; Leergaard et al. 2012; Partzsch and Schüffny 2012; Supek and Magjarevic 2011) have demonstrated that biological networks presented small-world properties (Bassett and Bullmore 2006; Kaiser 2008; Sporns and Zwi 2004). A small-world network constitutes a compromise between random and nearest neighbor regimes resulting in a short average path length despite the predominance of local connections (Watts and Strogatz 1998), which reflects the high efficiency of the network in transmitting and processing information. However, there is still a lack of direct and systematic evidence to prove its effects on information encoding performance.

Traditional approaches to neural coding hypothesized that information was represented in the precise timing of individual spikes, the firing rate or inter-spike interval of individual neurons (Bialek et al. 1991; Shadlen and Newsome 1994). Subsequently, researchers found synchrony and oscillatory responses in neural population of visual cortex (Eckhorn et al. 1988; Gray and Singer 1989; Toyama et al. 1981) as well as in other sensory cortexes (Laurent and Naraghi 1994; Neuenschwander and Singer 1996; Pantev et al. 1991), and also revealed that the interaction and oscillation strengths from the primary visual cortex (V1) mostly depended on the outside stimuli (Samonds et al. 2006). However, those studies only focused on revealing the neural mechanisms through reporting experimental phenomenon deprived from neuronal data recorded from various animals, e.g., rats, cats, monkeys, etc.

In order to better understand the information encoding mechanisms of the brain, more and more researchers attempted to establish population encoding model (Haslinger et al. 2013; Pillow et al. 2008) or information integrity model (Wang et al. 2011) to describe the computing process of the local neuronal population. These studies demonstrate the important role of synchrony or oscillation in population encoding and information decoding. However, they did not take the connectivity structure of neuronal population into account. There are also some studies focusing on building models to explore the important role of the topological structure in representing the dynamics of neural population (Pernice et al. 2011, 2013; Trousdale et al. 2012, 2013). However, most of those models are established based on large quantity of neurons, which could not be examined with electrophysiological data.

In order to investigate the effect of small-world structure of neural population of V1 on population encoding and information decoding, we recorded neuronal population from rat V1 with microelectrode array (MEA) and analyzed the dynamics of the neuronal population under a small world-based framework (SWF for short) on the one hand. The electrophysiological results show that the variation between population responses to different topological shape stimuli characterized with SWF were significant than that not taking the small-world structure into account. On the other hand, we established a small world-based generalized linear model (SW-GLM for short) to describe the computation process of neuronal population in V1. The parameters were obtained from the electrophysiological data recorded from rat’s V1. The theoretical modeling results show that both the precision of predicted neuronal population activity and the decoding performance were improved by taking into account the small-world structure. Therefore, we provided both electrophysiological and theoretical evidence for demonstrating the effects of small-world structure on information encoding performance in the primary visual cortex.

Methods

Animal preparation

Our data were recorded from five Long Evans rats (2–4 months, 200–300 g) supported by Animal Center of Zhengzhou University. Each rat was anesthetized with urethane (1.3–1.5 g/kg) and restrained in a stereotaxic apparatus (David Kopf Instruments). The rectal temperature was monitored and maintained at 37.5 °C via a heating pad. The right eye was fixed with a metal ring to prevent eye movement. A craniotomy (diameter 3 mm) was performed over the left V1 and the dura was removed to permit the electrodes to be implanted. After this procedure, no new recordings were made for a period of at least 20 min.

Visual stimuli and electrophysiological recordings

Visual stimuli were generated with a PC containing a Matlab toolbox (Psychtoolbox) and displayed on an LCD (liquid crystal display) monitor (250 × 400 mm, refresh rate 60 Hz) positioned at a distance of 15 cm in front of the rat’s right eye. The receptive fields (RF) were first mapped with sparse noise stimuli, with single bright square (6.6° × 6.6°) flashing on a black background at each of the 12 × 12 positions in a pseudorandom sequence (20 flashes/position), applying reverse correlation (Jones et al. 1987). All the following stimuli were located covering the multiple RFs simultaneously. The orientation selectivity for each neurons were measured under drifting sinusoidal gratings with 12 different orientations repeated for 20 times (30° in steps, the time frequency of 2 Hz, the spatial frequency 0.1 cpd, the contrast 100 %). The synchrony strengths were calculated under this type of stimuli to evoke maximum response for most neurons. The simulated model was also examined under grating stimuli.

In order to analyze the response dynamics of recorded population, we designed another set of stimuli, i.e., figures with different shapes identity (the areas in white remain the same for all figures. The size and location of shape figures were adjusted to cover most RFs of the population, see Fig. 1. Each shape figure stimulus was randomly represented 50 times, 1 s in duration with a 1-s blank window interleaved between stimulus presentations to avoid adaptation.

Fig. 1
figure 1

Different shaped figures with the fitted RFs of the example dataset presenting behind. The right two are topologically equivalent and both are topologically different from the left one. The spatial profile of each recorded unit was indicated by gray ellipses shown behind, among which the one shown with thick ellipse was the example neuron referred to in Figs. 4a and 7a

Full size image

The neuronal data were all recorded with an MEA (4 × 4, intercontact distance 350 μm for adjacent electrodes). The signals were collected with a Cerebus system (Blackrock Microsystems) and amplified 4000 times and band-pass filtered between 250 Hz and 5 kHz. Both the timestamps and the waveforms of the detected spike were recorded (30,000 Hz sampling frequency) for applications of off-line data analysis.

Data analysis: small world-based framework to quantify population response

The population responses to various stimuli were quantified with a small world-based framework, whose schematic chart is shown in Fig. 2. The computing steps are summarized as follows:

Fig. 2
figure 2

The schematic description of the small world-based framework. The detailed description of computational process of SWF is given in “Methods”. The terms referred to in the figure are all consistent with the main text

Full size image

Firstly, spike sorting was performed to isolate single units applying an unsupervised clustering algorithm (Quiroga et al. 2004). In our experiments, we only took into account the neurons responding significantly to stimuli designed in our paper. To this end, we computed a baseline activity using the mean firing rate within 0.5–1.0 s after the stimulus onset, and defined a threshold as five times SD above the baseline. For each stimulus, the subthreshold (or spike) response was considered to be significant if the peak amplitude was above the threshold (Zhu and Yao 2013). The orientation tuning curve for each selected unit was drawn and its corresponding orientation sensitivity index (denoted by OSI) was calculated with the following equation (Niell and Stryker 2008):

$${\text{OSI}} = \left( {R_{\text{pref}} - R_{\text{org}} } \right)/\left( {R_{\text{pref}} + R_{\text{org}} } \right),$$
(1)

where R pref and R org represent the mean firing rate under gratings of the preferred and its orthogonal orientation, respectively. R org = (R pref−π/2 + R pref+π/2)/2.

Secondly, the response of each individual neuron was quantified in the following way (upper path in Fig. 2). The firing spike train of each candidate unit was time binned (10 ms/bin) and the post-stimulus time histogram (PSTH) was obtained by applying the spike trigged average (STA) method. The response index (denoted by RI) for each unit under each type of stimuli was calculated with

$${\text{RI}} = 10 \times \left| {\log_{10} (R_{|s|} /R_{|b|} )} \right|,$$
(2)

where R |s| and R |b| represent the mean firing rate within 0–0.4 s after stimulus onset and the firing rate of the last 0.4 s of the blank stimuli, respectively. A logarithmic transformation (Ringach et al. 2002) was performed to account for inhibitory units, whose R |s| was always less than R |b|.

Thirdly, the synchrony strengths between pairwise units were calculated with normalized cross-correlation histogram (CCH, statistically equivalent to shift predictor), which would remove the synchrony caused by the covariations of the firing rates resulted by the simultaneous presentation of the stimulus and retain the synchrony caused by the intrinsic interaction between pairwise cells (Aertsen et al. 1989). In our experiments, the peak value at ~0 ms time lag in the normalized CCH was considered as a synchrony of strength I ij between neurons i and j.

Fourthly, a small-world network was estimated based on the calculated synchrony matrix, preparing for quantifying the connection strength of the population (lower path in Fig. 2). The synchrony matrix was converted into a binary graph (called a connection matrix, denoted by W) by setting a proper threshold (denoted by Tr in Eq. 3).

$$\omega_{ij} = \left\{\begin{array}{ll} 1& \quad {I_{ij} \ge {\text{Tr}}} \\ 0& \quad {I_{ij} < {\text{Tr}}} \\ \end{array}\right.,$$
(3)

where 1 refers to connected and 0 unconnected. In our experiment, the threshold ensures to emerge a network with the largest small-worldness (see the next section).

Finally, the population response was quantified using Eq. 4, composed of two terms. The first indicates the response strengths contributed by individual units, while the latter term represents the coupling strengths of the population and relies mostly upon the connection structure, denoted by W (ω ij , i = 1,…, N; j = 1,…, N).

$$R_{\text{CA}} = \frac{1}{N}\left( {\sum\limits_{i = 1}^{N} {\alpha_{i} \left[ {\text{RI}} \right]_{i} + \beta \sum\limits_{i > j}^{N} {\sum\limits_{j = 1}^{N} {\omega_{ij} I_{ij} \left[ {\text{RI}} \right]_{i} \left[ {\text{RI}} \right]_{j} } } } } \right)$$
(4)
$$\beta = N/\sum\limits_{i > j}^{N} {\sum\limits_{j = 1}^{N} {\omega_{ij} } }$$
(5)

In Eq. 4, α i as the weighted value for cell i was determined by its orientation selectivity OSI (Eq. 1) (see “Results” for explanation). β is the contribution factor of connection strengths of the population and varies with the amount of connectivity existing in the population to make the output of models with different connectivity comparable. In our experiments, we introduced a control approach, which was based on an “all to all” fully connected network and was comparable to taking all possible connections among the population into consideration [denoted by full network-based framework (FNF)]. The parameters for each method are shown in Table 1.

Table 1 The parameters of the proposed framework compared with those of the control methods
Full size table

Data analysis: quantify the small-worldness of a network

We quantified the small-worldness of each thresholded network in the following way.

The average path lengths L (Eq. 6) and clustering coefficient C (Eq. 7) for the estimated network were first computed and then compared with those of 100 random networks, which preserved the same number of nodes and connections as the actual network and had nodes connected to each other randomly.

$$L = \frac{1}{{\frac{1}{2}N\left( {N - 1} \right)}}\sum\limits_{i > j} {d_{i,j} }$$
(6)
$$C = \frac{1}{N}\sum\limits_{i = 1}^{N} {\frac{{2E_{i} }}{{k_{i} \left( {k_{i} - 1} \right)}}} ,$$
(7)

where N is the number of nodes in the network. \(d_{i,j}\) indicates the distance between nodes i and j, which is defined as the number of nodes forming the shortest path connecting nodes i and j. E i represents the number of actual connections among node i and its neighbors. Let λ indicate the ratio between the average path length of the actual network L s and that of random networks L r , and γ the ratio between the average clustering coefficients of the two networks. That is, λ = L s /L r , γ = C s /C r . We computed the index S w by

$$S_{\text{w}} = \, \gamma /\lambda$$
(8)

A small-world property was inferred if λ ≈ 1 and γ > 1, or, equivalently, if the ratio S w > 1 (Achard et al. 2006). Therefore, we preferred the network with the largest Sw to represent the connection structure of the recorded neuronal population.

Small world-based generalized linear model

A small world-based generalized linear model was built here to provide theoretical evidence for demonstrating the effect of small-world network on encoding and decoding performance. The general schematic chart of the SW-GLM is shown in Fig. 3. The general idea of GLM model assumed that the spiking probability of each individual cell is related to the modulation by the stimulus, the past activity of the neurons including its own and the coupled ensemble’s (Okatan et al. 2005; Truccolo et al. 2005). Therefore, the spiking of individual neuron among the N-neuron group was always described with three sets of filters: a spatiotemporal stimulus filter (indicated with F); a post-spike-dependent filter (represented with P), which captures mostly the refractoriness, recovery periods and adaptation of the neuron; and a set of coupling filters C, which capture dependencies on the recent spiking of the coupled cells. The summed filter responses were then exponentiated to produce an instantaneous firing rate for each individual neuron. In our paper, we made a hypothesis that the small-world structure of the neuronal network would also affect the firing activities of the population. To this end, we fit the structure of the model to the recorded neuronal network.

Fig. 3
figure 3

Model schematic for the SW-GLM. Each neuron has a stimulus filter, a post-spike filter and coupling filters that capture dependencies on spiking in other neurons. And the coupled neurons constitute a small-world network, derived from the recorded neuronal population. Summed filters output passes through an exponential nonlinearity to produce the instantaneous spiking rate

Full size image

As is shown in Fig. 3, each neuron is modeled with three sets of filters with those coupled neurons constituting a small-world network. Given the stimuli s, and the spiking history \(H_{1:i}^{1:N}\) of an N-neuron population, the conditional intensity function λ of a single cell’s spiking activity in the ith time tin t i is

$$\lambda \left( {t_{i} \left| {H_{1:i}^{1:N} ,s,\theta } \right.} \right) = \lambda_{F} \left( {t_{i} \left| {s,\theta_{F} } \right.} \right)\lambda_{P} \left( {t_{i} \left| {H_{1:i} ,\theta_{P} } \right.} \right)\lambda_{C} \left( {t_{i} \left| {H_{1:i}^{1:N - 1} ,\theta_{C} } \right.} \right),$$
(9)

where \(\theta = \left\{ {\theta_{F} ,\theta_{P} ,\theta_{C} } \right\}\) indicates the parameters of each filter to be estimated from the physiological data, \(\lambda_{F} \left( {t_{i} \left| {s,\theta_{F} } \right.} \right)\) is the spiking intensity induced by the extrinsic covariate of the stimuli s and is modeled as \(\lambda_{F} = \exp \left( {F \cdot s} \right)\). The corresponding filter F was represented with spatial and temporal filters:

$$F\left( {x,y,\theta ,\tau } \right) = F_{\text{s}} \left( {x,y,\theta } \right)F_{\text{t}} \left( \tau \right)$$
(10)

with \(F_{\text{s}} \left( {x,y,\theta } \right)\) denoting a spatial filter and obtained by approximating the measured RF with 2D Gabor function (Eq. 11) in a least-squares way, and F t (τ) the temporal filter determined by fitting the temporal function of RF response amplitude with gamma distribution function (Eq. 12).

$$F_{\text{s}} \left( {x,y,\theta } \right) = Ae^{{ - \left( {(x - x_{0} )\cos \theta + (y - y_{0} )\sin \theta } \right)^{2} /\sigma_{x}^{2} - \left( {(y - y_{0} )\cos \theta - (x - x_{0} )\sin \theta } \right)^{2} /\sigma_{y}^{2} }}$$
(11)
$$F_{\text{t}} \left( \tau \right) = U\left( {\tau - \tau_{0} } \right)^{\alpha } e^{{ - \left( {\tau - \tau_{0} } \right)/\sigma + A_{0} }} ,$$
(12)

where A is amplitude of the RF response, \(\sigma_{x}\) and \(\sigma_{y}\) determine the width and length of the RF, θ is orientation of the RF, at each temporal delay. In Eq. 12, \(U,\alpha ,\sigma ,\tau_{0}\) and \(A_{0}\) are free parameters.

\(\lambda_{P} \left( {t_{i} \left| {H_{1:i} ,\theta_{P} } \right.} \right)\) and \(\lambda_{C} \left( {t_{i} \left| {H_{1:i}^{1:N - 1} ,\theta_{C} } \right.} \right)\) of Eq. (9) represent the components of the intensity function conditioned on the spiking history \(H_{1:i}\) of its own and that of other K − 1 coupling neurons. They are, respectively, modeled as:

$$\lambda_{P} = \exp \left[ {\mu_{0} + \sum\limits_{n = 1}^{Q} {p(t_{n} )H_{i - n} } } \right]$$
(13)
$$\lambda_{C} = \exp \left( {\sum\limits_{k = 1}^{N - 1} {\omega_{k \cdot } \cdot c_{k} \cdot H^{k} } } \right) = \exp \left( {\sum\limits_{k = 1}^{N - 1} {\omega_{k \cdot } \cdot \sum\limits_{n = 1}^{{Q^{'} }} {c_{k} \left( {t_{n} } \right)\,H_{i - n}^{k} } } } \right).$$
(14)

In Eq. (13), \(\mu_{0}\) is the cell’s baseline log-firing rate, \(p(t_{n} )\) indicates the gain coefficients of post-spike filter P at time \(t_{n}\). In Eq. (14), \(\omega_{k \cdot }\) is the element of connection matrix (Eq. 3) and indicates whether the target neuron was connected with the kth neuron. \(c_{k}\) and \(H^{k}\) represent the coupling filter and the history spiking of the kth coupling neuron, respectively. The post-spike filter P and coupling filter C were represented using a set of raised cosine basis of the form

$$b_{j} \left(t \right) = \left\{\begin{array}{l} {\frac{{\cos \left( {a\log \left[ {t + d} \right] - \phi_{j} } \right) + 1}}{2}}\;\;\;{{\text{If}}\,\,a\,\log (t + c) \in \left[ {\phi_{j} - \pi ,\phi_{j} + \pi } \right]} \\ 0 \qquad\qquad\qquad\;\;\quad{\text{Otherwise}}\\ \end{array}\right.$$
(15)

All model coefficients \(\theta\) and their significance p values were estimated with physiological data using standard maximum-likelihood methods (Paninski 2004). To make a control analysis, we designed a control model without considering the effect of the network structure (called full-GLM), which means the third part was modeled as:

$$\lambda_{C}^{\text{full}} = \exp \left( {\sum\limits_{k = 1}^{N - 1} {c_{k} \cdot H_{k} } } \right) = \exp \left( {\sum\limits_{k = 1}^{N - 1} {\sum\limits_{n = 1}^{{Q^{'} }} {c_{k} \left( {t_{n} } \right)H_{i - n}^{k} } } } \right).$$
(16)

The parameters in Eq. 16 have the same meaning with in Eq. 14.

Verifying the encoding effect of the SW-GLM

The model was validated using a novel 1-s stimulus repeated for 20 presentations. The raster was obtained and the time-varying average response (PSTH) was computed in 1-ms time bins, smoothed with a Gaussian kernel of widths 2 ms. In order to quantify the prediction accuracy of the model, the residual error [defined as the difference between the true and predicted values (Andersen et al. 1992)] was computed and averaged over nonoverlapping moving time windows [t T−I  − t I ] (\(T - I \ge 1\)) using Eq. (17).

$$e\left( {t_{T} } \right) = \left| {\sum\limits_{i = T - I}^{T} {B_{i} } - \int_{{t_{T - I} }}^{{t_{I} }} {\lambda \left( t \right){\text{d}}t} } \right|/M,$$
(17)

where B i denotes the mean firing rate of the neuron in the ith time bins, \(\lambda\) is the conditional density of the model indicating the predicted firing rate in each time bins, and M represents the number of the time bins.

In addition, we computed the Pearson correlation coefficient of the model’s responses between each pair of repeats and averaged the correlation coefficients over all pair-wise combinations to quantify the predicting reliability of the model.

Checking the decoding performance of SW-GLM

We applied the regularized logistic regression (Bishop 2007) to decode the stimulus orientation from single-trial population response using 240 different orientation grating stimuli. The decoder was to decide whether this population response occurred in a trial under the stimulus with orientation φ 1 or φ 2. Such a classification has proven useful to assess the quality of different population coding schemes (Berens et al. 2011). We trained the logistic regression model using the glmnet toolbox (Friedman et al. 2010) in MATLAB with L 1/2 regulation instead of L 1 regulation. The L 1/2 regulation has proved to achieve better constringency performance than L 1 regulation (Zhao et al. 2012). For each pair-wise combination of stimulus orientation, the cross-validation was performed using the repeated random subsample technique (80 % training data, 20 % test data) for the whole regularization path (50 regularization parameters spaced between e −10 and e). The correct percentage over all tested data for each combination of orientations was averaged to estimate the decoding performance based on population response predicted by the model.

Results

A total of five datasets recorded from five rats’ V1 were used in the paper. The detailed results of one example dataset are presented together with the key results of others.

Response dynamics of recorded neuronal population

We analyzed the dynamics of population response to visual stimuli by quantifying response of recorded neuronal population applying SWF and made a comparison with results of a control approach by taking all possible connections amongst the neuronal population into account (FNF).

The individual response was first analyzed. The Gaussian-fitted PSTHs of each isolated unit under gratings were mapped (Fig. 4a). From Fig. 4a we can see that the cell responded significantly after the stimulus onset and that the response lasted for hundreds of milliseconds. The orientation tuning curves of the strongly responding neurons were mapped (Fig. 4b) and the corresponding OSI for each neuron was then measured (Eq. 1). The individual response was computed using Eq. 2, and compared with the calculated OSI to see whether the orientation selectivity would contribute to the response strength for each individual neuron. Figure 4c shows the comparison results for all recorded neurons (n total = 72). The raw data were fitted with a linear function. From the comparison we can see that the two groups of data were positively correlated with each other, suggesting that the value of OSI could approximately represent the response potential to varied visual stimuli, and the contribution made by those more active units should be larger. Thus, to quantify the population response (Eq. 4), we just simply took OSI value as the weighted value α for each individual neuron.

Fig. 4
figure 4

Response dynamics of individual neurons. a The raster (upper) and PSTH (below) of an example neuron. The solid red line shows the averaged firing rate and the shadow indicates the std error. The stimulus is repeated 20 times. The dashed line marks the timestamps of stimulus occurrence. b Orientation tuning curve for the example neuron. OSI = 0.5880. c The response indices (represented with RI) of each unit from the 5 datasets are plotted against their orientation selectivity indices (indicated with ORI). The best linear fit (solid line) is shown in addition to the goodness-of-fit statistics (denoted with r) and the p value

Full size image

The synchrony strengths between pair-wise neurons were then calculated with the normalized CCH, derived from the normalized joint-PSTH, to construct the synchrony matrix (Fig. 5a). A proper threshold was chose to produce a network with small-worldness (denoted by S w) as large as possible. The thresholded binary connection matrix (S w = 1.4354) is shown in Fig. 5b.

Fig. 5
figure 5

The synchrony matrix (left) and connection matrix (right) of the recorded neuronal population

Full size image

The population responses to visual stimuli were finally quantified using Eq. 4. To examine whether the characterized population response was an efficient way to discriminate visual stimuli, we performed a post hoc pairwise comparison between groups of response to different shape stimuli. The horizontal blue line indicates the distribution of population response (Fig. 6a) for 50 repetitions of “cross” stimulus (Mode 1 shown in Fig. 1). Blue line denotes in the target group, prepared to be compared with the other 2 groups of responses. If any of them had mean significantly (p < 0.05) different from the target group they would be present in red (as in Fig. 6a), otherwise shown in gray (as in Fig. 6b). To further quantify the variation between different groups of data, we performed one-way analysis of variance statistics (significance p values with a 95 % confidence interval) and calculated effect size [Hedge’s g (Fritz et al. 2012)]. The results are summarized in Table 2. From both Fig. 6 and Table 2, we can see that the population responses to similar topological shapes do not present any separation (Mode 2 versus Mode 3) for both methods (SWF and FNF). But the observed difference is significant in SWF-based population response to different topological shapes (Mode 1 versus Mode 2, and Mode 1 versus Mode 3) (p < 0.001 and Hedge’s g > 0.8) and not significant (p > 0.05 and Hedge’s g < 0.5) for FNF, the control methods without considering the small-world network structure.

Fig. 6
figure 6

A post hoc pairwise comparison between population responses to different shape identity. a The population responses were calculated under SWF. b The population responses were calculated under FNF

Full size image
Table 2 The comparative results of one-way analysis of variance (ANOVA) test and the calculated effect size
Full size table

Parameters of the SW-GLM

As mentioned in “Methods”, the established model mainly composed of three sets of filters.

First, the spatiotemporal RF was first mapped for each neuron. The spatial profile and the temporal dynamic of the RF were determined using Eqs. 11 and 12 to estimate the spatial and temporal stimulus filters. Figure 7a shows the measured RF of an example neuron at different time delays (upper) and the response amplitude as a function of time (lower, indicated with black empty dots), which was fitted with gamma distribution function to construct the temporal stimulus filter. The estimated spatial filters of the 15-neuron group (Fig. 1) formed an approximate complete mosaic covering a small region of visual space.

Fig. 7
figure 7

The fitted filters for the SW-GLM. a Illustration of parameter fitting of spatiotemporal RF. Spatial RF in each time interval (above each image) after the onset of the stimuli was fitted with 2-D Gabor function (upper, white ellipse). Shown are spatial RF profiles at nine temporal delays (black circles in lower trace). The second frame shown is the RF at the time of peak response. Amplitude of Gaussian fit is plotted against time delays (lower trace) to obtain the temporal stimulus filter. b Seven-dimensional basis for post-spike filters (upper) and four-dimensional basis for coupling filters (lower). c The estimated post-spike filter. d The coupling filters from the coupled neurons with the example neuron

Full size image

Second, the post-spike-dependent filter and the coupling filter were both represented using a set of raised cosine basis (Eq. 15, Fig. 7b) according to the inter-spike interval distribution and the temporal structure of the normalized cross-correlation histogram. The temporal range was chosen after the observation that the magnitude of most estimated filters decline back to zero well within the first 40–60 ms. The exponentiated post-spike and coupling filters are shown in Fig. 7c, d.

The predicted response of the established model

The predicted accuracy was first compared between SW-GLM and full-GLM to assess the effect of the network connectivity on neural encoding. The response to 20 repeats of novel stimuli was predicted under the SW-GLM, and the mean firing rate in 1-ms time bins as a function of time (PSTH) was smoothed with a Gaussian kernel of widths 2 ms. The residual error between the predicted PSTH and that of true recorded neuron is shown in Fig. 8a, compared with the predicted residual error of the control model. In order to quantify the predicting effect, we calculated the predicted accuracy of each model using Eq. (17). The comparing results across all recorded neurons (n total = 72) from five datasets are shown in Fig. 8b, from which we can see that SW-GLM predicted more accurate response (37.6 ± 6.2 % less residual error) than the full-GLM. The comparing result suggests that the computation described by the SW-GLM model is more suitable to describe encoding process of the true neurons recorded in V1.

Fig. 8
figure 8

The examination results under SW-GLM compared with the full-GLM. a The PSTH of the target neuron to 1-s novel stimulus in 1-ms time bin (raw PSTH), the predicted response under SW-GLM, comparing with that under full-GLM. b The predicted error of the full-GLM for all 72 neurons plotted against those of the SW-GLM. The one plotted with solid black circle was the predicted result for the target neuron. c The estimated reliability of response predicted by the full-GLM for all 72 neurons plotted against those under the SW-GLM. The SW-GLM predicts more reliable responses than the control model. The one plotted with solid black circle was the predicted result for the target neuron. d Average decoding performance for both numerical models plotted as a function of Δφ

Full size image

We then measured and compared the reliability of the response predicted by both models. Across all 72 neurons recorded from five rats’ V1, SW-GLM predicts much more reliable response (p < 10−4, Wilcoxon signed-rank test, Fig. 8c) than the control model. The comparing result indicates that the small-world structure may help to produce less noisy spikes.

The decoding performance of SW-GLM

However, the above comparisons only suggest that the small-world structure of the neuronal population might be necessary to predict responses of V1 neurons, but not enough to extract the stimulus information that the responses carry. To address this issue, we used the L 1/2 regularized logistic regression to decode the grating orientation from the response predicted by the SW-GLM and compared the decoding performance with that of the control model. The percent correct between the 100–500 ms after stimulus onset were averaged to present the decoding accuracy for each combination of grating orientation. Figure 8d shows the comparison result of the decoding accuracy plotted against \(\Delta \varphi\) (\(\Delta \varphi = \left| {\varphi_{1} - \varphi_{2} } \right|\)). The comparison result indicated the decoding performance was largely improved by taking the topological structure of the neuronal population into account, which suggested that the network structure that neurons gather together is critical for extracting the sensory information contained in the population response.

Discussion

We investigated the role of the small-world network that emerged from the neuronal population in V1 in encoding and decoding outside information with electrophysiological and numerical simulating experiments. Both aspects of results demonstrated our hypothesis that the small-world network structure plays an important role in the process of population encoding and information decoding, which is consistent with the conclusions that computations performed by visual cortex involve the activities of individual single units as well as the network connectivity (Bock et al. 2011).

Physiologically, we quantified the response of neuronal population recorded from rat’s V1 applying a small world-based framework and made a comparison with the control approach without considering the small-world structure emerged from the neuronal population in V1. The small-world property has already been proved to exist in anesthetic cat V1 by Yu et al. (2008) and in alert monkey V1 by Gerhard et al. (2011). The comparison results show that the SWF-based responses were more efficient to discriminate different shape identities compared to those by simply summing up all synchrony strengths between pairwise neurons. Coincidentally, shape processing in the visual cortex has been reported to start from V1 along the ventral visual stream (Ungerleider and Mishkin 1982) and contour perception has been linked to lateral interactions between neighboring neurons in V1 (Loffler 2008). The control approach was prevalently used to analyze characteristics of population response. For example, Zhou et al. (2008) draw conclusions that the population response would be modulated with the spatial integrity of the visual stimuli by summing up synchronous strengths between each pair of neurons. In conclusion, the physiological results not only verify the effect of the small-world network structure on encoding global visual shape information, but also provide us a general tool for analyzing neuronal data of much larger quantity, such as two-photon imaging data (Matsui and Ohki 2013), which may help to better understand neural circuit and its development.

Mathematically, we simulated the computation of neuronal population with a small world-based GLM model. Both the predicted accuracy and decoding performance were improved compared to the control model without a small-world structure, suggesting the significant role of the small-world structure in the process of encoding and decoding information. Similarly, Pillow et al. (2008) have used a GLM model to prove the role of correlated activity in the retinal coding of visual stimuli, but they did not consider the influence of network structure. Pernice et al. (2013) and Trousdale et al. (2013) have, respectively, established numerical models and reported the strong influence of the network structure on the population activity dynamics, while their models were also based on large of quantity of neurons (more than 1000), which out of the recording ability of current available MEA. In addition, Zheng et al. (2014) have established a small-world neuronal network model with each neuron described with a two-dimensional Rulkov map neuron and proved the introduction of small world connectivity would enhance the temporal order of the spiking activity.

The study carried here suggests the small-world structure may convey additional information about the sensory stimuli with simulation and experiment, which provides a general framework for understanding the importance of structured synchrony activity in populations of neurons. However, the connectivity obtained from our experimental data or simulating model may not be strictly consistent with the real intrinsic connectivity of visual cortex. To further match it with the visual cortex’s anatomical connectivity would need some new neuroinformatic tools and some systematic multi-scale approaches, which will be one of the tasks in our future research.