OPENSUBJECT AREAS:
NEURAL ENCODING
NEUROINFORMATICSNETWORK MODELS

Simultaneous silence organizes
structured higher-order interactions inneural populations
Hideaki Shimazaki1, Kolia Sadeghi2, Tomoe Ishikawa3, Yuji Ikegaya3,4 & Taro Toyoizumi1,5

NEURAL CIRCUITS
1

Received
13 December 2014
Accepted
18 March 2015
Published
28 April 2015

Correspondence andrequests for materials
should be addressed to
H.S. (shimazaki@
brain.riken.jp) or T.T.
(taro.toyoizumi@
brain.riken.jp)

RIKEN Brain Science Institute, 2-1 Hirosawa, Wako-shi, Saitama 351-0198, Japan, 2Commonwealth Computer Research Inc.,
1422 Sachem Pl., Unit #1, Charlottesville, VA 22901, U.S.A., 3Graduate School of Pharmaceutical Sciences, The University ofTokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan, 4Heart to Heart Science, Center for Information and Neural Networks, 1-4
Yamadaoka, Suita, Osaka, 565-0871, Japan, 5Department of Computational Intelligence and Systems Science, Tokyo Institute ofTechnology, Yokohama 226-8502, Japan.

Activity patterns of neural population are constrained by underlying biological mechanisms. These patterns
are characterized not only by individual activity rates and pairwise correlations but also by statisticaldependencies among groups of neurons larger than two, known as higher-order interactions (HOIs). While
HOIs are ubiquitous in neural activity, primary characteristics of HOIs remain unknown. Here, we reportthat simultaneous silence (SS) of neurons concisely summarizes neural HOIs. Spontaneously active neurons
in cultured hippocampal slices express SS that is more frequent than predicted by their individual activity
rates and pairwise correlations. The SS explains structured HOIs seen in the data, namely, alternating signs
at successive interaction orders. Inhibitory neurons are necessary to maintain significant SS. The structured
HOIs predicted by SS were observed in a simple neural population model characterized by spiking
nonlinearity and correlated input. These results suggest that SS is a ubiquitous feature of HOIs that
constrain neural activity patterns and can influence information processing.

I

nformation in the brain is represented by the collective spiking activity of multiple neurons1. Activity patternsof observed neurons are highly structured due to various underlying biological mechanisms including direct
anatomical connections2,3, indirect connections mediated by unobserved neurons4,5, and intrinsic nonlinearity
of individual neurons6,7. However, exploration of this structure is non-trivial due to limited data size in comparison to possible combinations of activity patterns that grow exponentially with population size.
To infer the structure of neural activity patterns from limited amount of data, the maximum entropy principle
has been successfully applied8,9. Under this principle, the probability distribution of activity patterns is estimatedto be the least structured distribution that is consistent with a set of observed activity statistics. Conventionally
this maximum entropy distribution is statistically characterized by parameters of different orders, where theorders refer to the numbers of subset neurons that these parameters constrain. The model with the first-orderparameters fits to the observed activity rates of individual neurons. The model that additionally includes thesecond-order parameters further adjusts the observed deviations of pairwise correlations from the chance coincidence expected from the individual activity rates. The second-order parameters are referred to as pairwiseinteractions. More generally, the model that includes up to the k-th (k 5 3, 4, . . . ) order interactions adjusts
simultaneous activation rates of k neurons from the expectation based on interactions up to the (k-1)-th order.
Interactions beyond the pairwise interactions (k. 2) are collectively termed higher-order interactions (HOIs)10,11.
Notably, these interactions refer to statistical dependency of neurons, and do not necessarily involve anatomical
connections.
In earlier studies, individual activity rates and pairwise correlations alone could explain , 90% of variability inactivity patterns of small populations of retinal ganglion cells8,9 and cortical neurons12,13. However, this does not
exclude the existence of HOIs or limit their contribution to information processing. Indeed, the addition of HOIs
to a statistical model significantly improved the goodness-of-fit to neural activities obtained from multi unit
activity14,15, single unit activity5,16–20, and local field potential21,22 in both in vivo and in vitro preparations.
Furthermore, HOIs are relevant in neural information coding14,16,18,23. However, previous studies have not
identified a key feature in HOIs that summarizes the principal role of seemingly diverse HOIs.

SCIENTIFIC REPORTS | 5 : 9821 | DOI: 10.1038/srep09821

1www.nature.com/scientificreports
One of the most striking features of neural population activity issimultaneous silence (SS). The spiking activity of individual neurons
is known to be sparse24. As a result, the most commonly observed
activity pattern in typical networks is the pattern in which all neurons
are silent. Does SS involve HOIs? Indeed, departures from the level of
expected SS from individual activity rates and pairwise correlations(excess SS) were empirically reported previously16–20. However, thesignificance of SS in characterizing HOIs of the population activity is
not well understood.
Here, we examine SS in population activity of the hippocampal
CA3 networks in cultured slices. Previous studies demonstrated
that CA3 pyramidal cells in the organotypic slice cultures are wiredwith an in vivo-like connection probability of 15–30%3, and theirspontaneous spike rates are closer to those of in vivo hippocampal
neurons25, compared to neurons in acute slice preparations. We
demonstrate that most local groups of hippocampal neurons thatpossess HOIs express excess SS. A single parameter that quantifies
SS accounts for about 20% of the variability in population activity
patterns that is produced by numerous HOIs. We then confirm
specific oscillatory structure of HOIs at successive interaction
orders predicted from the SS. Through modeling, we also demonstrate that correlated population activity caused by spiking nonlinearity and correlated input exhibits the same structure of HOIs,
and that this structure conveys information of input. These results
suggest that neurons are operating in a unique regime where they
are constrained to be silent simultaneously.

Results
Simultaneous silence and HOIs of hippocampal neurons. We
analyzed the spontaneous spiking activity of putative neurons in
the hippocampal CA3 area of organotypic slice cultures, measured
by the Calcium imaging method. Slices were prepared from postnatal
day 7, and then cultivated from day 7 to 14 (see Methods). Neuronal
activity was detected by onsets of calcium transients3,26–28, which
provided event-timing data with a resolution of 100 ms. Fig. 1A
and B display an example of population event activity of a single
slice culture, and spatial positions and activity rates of individual
neurons. We analyzed n~20 slices in total, and found thefollowing features. First, activity rates of the neurons from all slices
were distributed close to a log-normal distribution (Fig. 1C),
similarly to spike rates of in vivo hippocampal CA3 neurons of
awake rodents29,30. The rates of calcium events in individual cells
computed from 2122 neurons in 20 slices were 0.073 6 0.097(mean 6 standard deviation (SD) events/s; median 0.035,
interquartile range 0.01–0.097 events/s). Notably, activity rates ofneurons in cultured slices were close to those under an awake invivo condition25. Second, the activity of pairs of neurons was only
weakly correlated (Fig. 1D). Average correlation coefficient was
0.033 6 0.065 SD (detection in a 100 ms window). A crosscorrelogram revealed that, on average, the activity of pairs ofneurons was not correlated after a , 400 ms timelapse (Fig. 1D
inset). Third, intracellular voltage recordings under the same
experimental conditions all reveal uni-modal distributions ofmembrane potentials (Fig. 1E). Hence, no obvious sign of asuperposition of UP and DOWN states was detected.
To analyze the correlated activity of multiple neurons, 50 groupsof N~10 neighboring neurons were selected from each of 20 slices(see an example group of neurons shaded in pink in Fig. 1B andevents marked in red in Fig. 1A), for a total of 1000 groups of 10
nearest-neighbor cells. The centers of groups were sampled according to the spatial density of cells in the CA3 area (See Methods). The
average ‘radius’ of the 1000 groups was 36.6 (613.4 SD) mm, where
the radius of a group was computed as the mean Euclidean distanceof its cell positions from the group’s center position. We then represented the activity of the ith neuron ði ~ 1,    ,10Þ in a time window
by a binary variable xi ~{0,1}, where ‘1’ denotes an active state in
SCIENTIFIC REPORTS | 5 : 9821 | DOI: 10.1038/srep09821

which at least one event occurred, and ‘0’ represents an inactive, or‘silent’, state in which no events occurred (Fig. 1F). We used a 400 ms
time-window in the subsequent analyses to incorporate the temporalcorrelation observed in the cross-correlogram (c.f. the shaded interval in Fig. 1D inset).
To examine if hippocampal neurons exhibit collective activity
beyond what can be explained by pairwise interactions, we compared the activity patterns of a group of observed neurons with
8,9,31,
those predicted frommaximumhX a pairwise
i entropy modelX:::
hi x i z
hij xi xj . This model provides
p2 (x1 , ,xN )* expiivj

the least structured probability distribution that is consistent
with the observed activity rates of individual neurons and correlations between pairs of neurons. The parameters fhi ,hi jg were
adjusted to fit these statistics. We call this model a pairwise
model hereafter. First, we examined if the neurons exhibited
SS beyond that predicted by the pairwise correlations. To thisend, we compared the observed probability of the pattern inwhich all of 10 neurons are simultaneously silent with its
probability according to the pairwise model. Figure 2A displays
a distribution of percentage deviation of observed SS
probabilities from the prediction of the pairwise model,
pdata ð0,    ,0Þ{p2 ð0,    ,0Þ, where pdata ð0,    ,0Þ is the observed
p2 ð0,    ,0Þ
probability of SS. In some groups, the pairwise model tended tounderestimate the occurrence probability of SS of 10 neurons.
This discrepancy has to be explained by HOIs in the data.
To examine the contribution of HOIs to population activity, we
computed the fraction of entropy that is explained by HOIs. This
fraction, referred to as the percentage entropy margin for HOIs, is
H2 {Hdata, where H2 is the entropy of thequantified as DHOI ~H2
pairwise model and Hdata is the entropy of the observed histogramof population activity patterns. We call Hdata the data entropy in thefollowing. The data entropy is characterized by all of the first, second,and HOIs. Therefore, the difference between H2 and Hdata must be
explained by HOIs. We found that the distribution of DHOI exhibiteda long tail (Fig. 2B). This indicates that there were a noticeable
number of groups in which HOIs played a much stronger role inshaping population activity. Finally, we explored the relation
between the contributions of HOIs to the probability of SS. We foundthat the groups expressing higher/lower probabilities of SS than the
pairwise model coincided with the groups possessing large entropy
margins for HOIs (Fig. 2C). The positive correlation between thesetwo values in Fig. 2C (Spearman’s rank correlation coefficient 0.69,
p ,0.001) implies that a significant portion of the HOIs of the CA3
neurons may be explained by the SS. The rank correlation coefficientwas higher (0.92) and statistically significant (p , 0.001) if we
analyze non-overlapping groups.
Simultaneous silence is a ubiquitous feature of HOIs. To directly
examine the contribution of the SS to the entropy explained by HOIs,
we constructed a maximum entropy model that augments the pairwise
model with a single additional term to account for the probability of SS
observed in the data. We refer to this model as the SS model:
Pss ðx1 ,    ,xN Þ*exp

hXi

hi xi zXivj

iN
hij xi xj zh0 Pi~1 ð1{xi Þ : ð1Þ

Here a single parameter, h0 , was introduced to account for the
probability of SS of N neurons. Positive or negative h0 indicates thatthe probability of SS of all neurons is more or less than predicted by the
pairwise model, respectively. Importantly, this new SS term isequivalent to adding specific structured HOIs into the pairwise
model. By expanding the SS term into the standard HOIXX
ð1{xi Þ~h0 {h0x zh0
coordinates, we obtain h0
i iivj

xi xj {h0PPx x x zh0
ivjvk i j k

iPxxx x:
ivjvkvl i j k lHence, increasing
2

0

200

400

600

800

1000

600
Y position [μm]B
Neurons

A

0.3
0.2

300

0.1
0

1200

0
0

Time [s]

D

0.7
0.6
0.5
0.4
0.3
0.2
0.1

Cross correlogram
2

20

1

0

10
-1

−4

10

−3

−2

10

−1

10

10

0

10

1

0.5

1

0

F

1.2

0.2

0.4

0.6

0.8

Correlation coefficient

Event sequences

Binary patterns

Bin size 400 ms

1000100000000
000100010000
000100000010
000100000000
000000000000
000100000000
000000010000
000000000000
000100000000
000100000010

1

0.8

Neurons

Probability density

0Time [s]Activity rate [Hz]

E

-0.5

0

0
10

600
X position [μm]30Probability density

Probability density

C

Activity rate [Hz]www.nature.com/scientificreports

0.6
0.4
0.2
0
−70

10
−65

−60

−55

−50

Membrane potential [mV]

686

688

690

Time [s]Figure 1 | Ensemble activity of CA3 putative neurons detected by Calcium imaging. (A) Ensemble activity of 45 neurons from a single hippocampal
slice. Small vertical ticks indicate events detected from calcium imaging signals. Ensemble activity of an example group of 10 neurons is marked in red. (B)Spatial distribution of neurons in the CA3 area of the slice in A. Each filled circle represents a position of a neuron. The color indicates activity rate of eachneuron. The pink area corresponds to the example group highlighted in A. (C) Distribution of activity rates from neurons in all 20 slices. Solid line is afitted log-normal distribution. (D) Distribution of correlation coefficients calculated from the event sequences (within a 100 ms window) from all thepairs of neurons in 1000 neighboring groups from 20 slices. The inset shows an average cross-correlogram from all the pairs of neurons. Dashed lines
indicate 6 2 SD of the correlogram at 1–2 sec lags. The gray shading (20.2 ms to 1 0.2 ms) indicates the interval where the correlgoram exceeded thedashed lines. (E) Distributions of membrane potentials recorded from neurons (n 5 7) in hippocampal slice cultures under the same condition as
described in Methods. Different colors indicate different neurons. In all cases, the densities of the membrane potentials were characterized by a unimodal
profile. (F) Construction of binary patterns from event sequences. The event sequences are binned using a window of 400 ms. In each bin, we denote ‘0’ ifthere is no event, and ‘1’ if there is at least one event.

or decreasing the total period of quiescence is equivalent tointroducing a single parameter to the HOIs with alternating signs
for different orders of interaction. In addition to capturing
individual activity rates and pairwise correlations, the SS model
explores this 1-dimensional structure in the high-dimensional spaceof HOIs to fit the rate of SS. We fitted the SS model to the same 1000
groups of 10 hippocampal neurons obtained from 20 slices. Note that,in each group, the fitted first and second order parameters of the SS
model are generally different from those of the pairwise model becauseof the newly introduced SS term.
We compared goodness-of-fit of the SS model with that of the
pairwise model (Fig. 3A). The ordinate of the panels represents
percentage differences between observed and predicted SS probabilities of sub-groups of r ð~ 1,    ,10Þ neurons by the two models (Left, the pairwise model; Right, the SS model). By definition,
the pairwise model adjusts the silence rates of individual neurons(equivalent to 1 minus activity rates, r~1) and pairs (r~2)(Fig. 3A Left panel). However, the pairwise model fitted to the
SCIENTIFIC REPORTS | 5 : 9821 | DOI: 10.1038/srep09821

data underestimated probabilities of SS for larger sub-groups ofneurons. This means that many sub-groups of hippocampal neurons expressed SS more often than chance as predicted from theiractivity rates and the pairwise correlations. In contrast, the SS
model additionally accounts for the probability of SS of all 10
neurons in a group (see the complete match of the data and prediction at r~10 in addition to r~1,2 in Fig. 3A Right panel).
Order of magnitude reductions in the differences were observedin the SS of many sub-groups (r ~ 3,    ,9). (Note the scale difference in the Left and Right panels.).
We tested the excess or paucity of SS using the SS model against anull hypothesis of no such activity (i.e., the hypothesis that the pairwise model is sufficient to characterize the data). Here, we used
x2 -tests11 with multiple comparison correction using theBenjamini-Hochberg-Yekutieli method with a false discovery rateof 0.05 to assess if the SS term significantly improved the fitting ineach group (See Methods). Of 1000 groups, 156 groups (16%) from
10 slices rejected the null hypothesis (Fig. 3B). We call these groups3www.nature.com/scientificreports

A

140

Number of groups

120
100
80
60
40
20
0
−2

0

2

4

Percentage deviation of SS
probabilities observed from predicted [%]B
Number of groups

120

80

40

0
0

2

4

6

8

10

12

14

C

Percentage entropy margin for HOIs, ΔHOI [%]Percentage entropy margin for HOIs, ΔHOI [%]

10

10

10

1

0

−1

−4

10

−3

10

−2

10

−1

10

0

1

10

10

Absolute percentage deviation of SS
probabilities observed from predicted [%]Figure 2 | Sub-groups of 10 hippocampal neurons exhibit longer periods ofSS than predicted from pairwise interactions. (A) Distribution of thepercentage deviation of the observed probability of SS of 10 neurons from thepdata ð0,    ,0Þ{p2 ð0,    ,0Þ
prediction of the pairwise model,. Positive
p2 ð0,    ,0Þ
values indicate more frequent SS in the data than predicted by the pairwise
model. (B) Histogram of percentage entropy margins for HOIs computed as
H2 {Hdata. (C) Dependency of the percentage entropy margin for HOIs on
H2
the percentage deviation of the observed from predicted probabilities of SS.
The same color indicates groups selected from the same slice culture. 5
outliers were excluded from the plots.

that exhibit excess or paucity of SS the SS groups. Statistical properties of the SS as well as non-SS groups were summarized in Table 1.
Indeed, most of the groups (68%, 133 groups out of 197) thatexhibited relatively large margins of entropy for HOIs (DHOI . 3%)
were the SS groups. Note that, at this point, each of the SS groups
could have had either significantly positive or negative h0 . It turned
out that 154 out of the 156 SS groups exhibited significantly positiveh0 (Fig. 3B Right inset). Thus, virtually all the SS groups expressedsignificantly larger probability of SS than the corresponding pairwise
model. In these groups, the total number of bins for which all neurons were quiet was larger than expected by the corresponding pairwise model. In other words, activity was confined to a smaller
number of bins. Hence, we conclude that the population activity of
most groups exhibiting HOIs (DHOI . 3%) was significantly sparse intime. Note that the observed fraction of SS groups was robust to thenumber of groups sampled from each slice but typically increased
with the size of these groups (Fig. 3C).
Finally, we examined the relation between the percentage entropy
margin for HOIs, DHOI , and the percentage of the HOI entropy
explained by the SS. The latter entropy was computed as
H2 {Hss, where Hss is the entropy of the SS model. Figure 3D
Dss ~H2
displays scatter plots of these values for all groups. As predicted from
Fig. 2C, we observe significant positive correlation between DHOI and
Dss (Spearman’s rank correlation coefficient 0.52, p ,0.001). The
H2 {Hss. This ratio
dashed lines are isoclines of a constant ratio a~
H2 {Hdata
describes the fraction of entropy explained by the SS in the entropy
margin for HOIs. As expected, the SS groups (filled circles) typically
had large a. Figure 3E displays a distribution of a for the groups that
expressed large margin for HOIs (197 groups with DHOI . 3%). In
these groups, the single higher-order parameter of the SS explained
18.3% (interquartile range, 4.7–31%) of the entropy for HOIs (Fig. 3C).
Since we have only added a single parameter in the high-dimensional
space of HOIs, this result implies that the SS comprises one importantcharacteristic of the HOIs.
In order to assess biases that may be caused by limited samples inour data sets, we repeated our analysis using two alternative data sets
(Supplementary Fig. S1 online). First, we analyzed only one half ofthe data by taking every other bin of the original population activity
patterns for each slice. Second, we analyzed bootstrapped population
activity patterns, where the same number of patterns as the original
data were resampled with replacement in each slice. These two data
sets contain less variations of population activity patterns than theoriginal data. For the both data sets, the fraction of SS groups was
smaller than the 16% found in Fig. 3B. The fraction of the HOIs
explained by SS also decreased to less than a half of 18% found in
Fig. 3E. Because we did not overestimate these quantities after subsampling and resampling, it is unlikely that our original estimation(16% exhibits significant SS; 18% of HOIs is explained by SS)
overestimated the fractions expected from a larger number of samples. In sum, the analyses confirm significant SS in the data, andpredict the presence of the alternating signs of HOIs, a possibility
we directly test now.
Alternating signs of HOIs predicted by SS. If SS is a major feature of
the HOIs, we expect to find HOIs whose signs alternate depending on
the orders of interaction (c.f. the expansion of the SS term). In order
to directly examine the structure of HOIs, we consider a simple
maximum entropy model that includes a single global parameterfor each order of HOIs:
phHOI ðx1 ,    ,x10 Þ* exp
hX
i
X
X10
X
h
hxz
hi1 i2 xix1 ix2 zx   xik ,
i i i
i vi
i vvi i1
k~3 k
1

2

1

k

ð2Þ

where hk (k~3,4,    ,10) is a single parameter for the kth order
SCIENTIFIC REPORTS | 5 : 9821 | DOI: 10.1038/srep09821

4www.nature.com/scientificreports
5

0.5

4

0.4

3

0.3

2

0.2

1

0.1

0

0

−1

−0.1
−0.2

Pairwise model1 2 3 4 5 6 7 8 9 10

1 2 3 4 5 6 7 8 9 10Size of sub-groups, rSize of sub-groups, r

84%

50

15

Number of groups

16%

100

C

SS groups

Total 1000 samples

SS groups

10
5

−1

Other groups

0

1

2

Significant θ 0

3

0
0

2

4

6

8

10

12

25
20
15
10
5
0
4

14

α=0.14α=0.27α=0.01

10

10

0

−1

10 12 14

HOI groups (ΔHOI>3%)

α=0.40

10

8

E
1

α=1.0

1

6

Group size, N

Cumulative distribtion

Percentage entropy margin for HOIs, ΔHOI [%]Percentage entropy margin for HOIs, ΔHOI [%]D

Total number of groupsn = 100n = 200
n = 1000
non-overlapping

0.75
α = H2 - HSS
H2 - Hdata

0.5

H2
HSS

H2 - Hdata

B
Number of groups

SS model−0.3−3

H2 - HSS−2

Fraction of the SS groups [%]Percentage deviationof SS probabilities [%]

A

0.25
Hdata

0
0
−4

10

−2

10

10

20

40

60

80

Proportion of SS in HOIs, α [%]0

Percentage entropy margin for SS, ΔSS [%]Figure 3 | Significant SS is observed in the groups of 10 neurons exhibiting HOIs. (A) (Left) Comparison of the observed SS probabilities of sub-groupsof r neurons with predictions of the pairwise model. Abscissa, the size r of sub-groups. Ordinate, the percentage deviation of observed from predicted
average SS probabilities of sub-groups of neurons, where the normalization divides the difference by SS probability predicted from a pairwise model. The
comparison was performed for all possible sub-groups of r~ð1,2,    ,10Þ neurons in the 1000 groups of 10 neurons. Whiskers represent 1.5 times the
distance from 25th to 75th percentile. Dots are outliers. (Right) Comparison of the observed SS probabilities of the sub-groups with predictions of the SS
model. Note the difference in the scales of the ordinates in the Left and Right panels. (B) Distribution of percentage entropy margins for HOIs for the SS
groups (red, 16%), i.e. groups for which the SS model showed significantly better fit than the pairwise model, and the rest of the groups (gray, 84%). The
right inset shows distribution of the parameter h0 of the SS groups. The distribution was heavily biased toward positive h0 , indicating prevalent excess SS.
(C) Group size dependency of the number of SS groups. Solid lines with different colors indicate different numbers of groups selected from each slice: 5,10, 50 groups per slice, for a total of n~100, 200, and 1000 groups of size N from 20 slices, respectively. The dashed black line is the result of selecting nonoverlapping groups from each slice (see Methods). Thus, the fraction of SS groups and its group size dependency were robust to the degree of overlap
between sampled groups in each slice. (D) Scatter plots of entropy margins for HOIs versus entropy margins for SS. The same color indicates groups
selected from the same slice culture. Filled circles indicate SS groups. Dashed lines represent different proportions of HOI entropy margin explained by theH2 {Hss
SS term, a~. We excluded 13 outliers from the plots. (E) Cumulative distribution function of the proportion of SS, a, in the groups exhibiting
H2 {Hdata
HOIs (DHOI . 3%).

HOIs. The term for the kth order interaction parameterized by aparameter hk sums all combinatorial interactions of k neurons
among 10 neurons. We call this model the homogeneous HOISCIENTIFIC REPORTS | 5 : 9821 | DOI: 10.1038/srep09821(hHOI) model. The hHOI model fitted to the data reproduces thehistogram of the number of active neurons in each time
bin14,19,20,32.
5

www.nature.com/scientificreports

Table 1 | Activity rates, correlation coefficients, and probabilities of SS computed from binary sequences using 400 ms window undercontrol and PTX conditions. Values are expressed as Mean (6 SD)

Control
SS groups
Non-SS gropus
PTX

Activity rate

Correlation coefficient

Prob. SS

0.039 (6 0.042)
0.032 (6 0.035)
0.040 (6 0.043)
0.027 (6 0.017)

0.060 (6 0.108)
0.172 (6 0.122)
0.040 (6 0.092)
0.920 (6 0.108)

0.728 (6 0.152)
0.831 (6 0.103)
0.706 (6 0.151)
0.965 (6 0.021)

In these data sets (DHOI . 3%), the hHOI model explained 24% of
the variability in population activity due to HOIs (interquartile range
H2 {HhHOI, where HhHOI is the entropy of
11–34%) as assessed by
H2 {Hdata
the hHOI model. This result indicates prevalent heterogeneity in the
HOIs. The result also upper bounds the fraction of entropy for HOIs
that could be explained by the single SS term. We then investigated
how much of this entropy is actually explained by the SS term.

β=0.2

10

1

B

β=0.4
β=0.6
β=0.8
β=1.0

HOI groups (ΔHOI>3%)

1

10

−1

−2

H2 - HSSβ = H2 - HhHOI

0.5

H2

0.25
HSS

H2 - HhHOI

10

0.75

H2 - HSS

10 0

Cumulative distribtion

Percentage entropy margin for hHOIs, ΔhHOI [%]A

Figure 4A displays relations between the percentage entropy margin
H2 {HhHOI, and by the SS
explained by the hHOI model, DhHOI ~H2
model, Dss . Similarly to Fig. 3D, the dashed lines are isoclines of
H2 {Hssb~, which quantifies the fraction of entropy explained
H2 {HhHOI
by the SS in the entropy margin for the homogenous HOIs.
Figure 4B shows a distribution of b for the groups exhibiting HOIs

HhHOI

0
0
10

−3

20

40

60

80

100

Proportion of SS in hHOIs, β [%]
−6

−4

10

10

−2

10

10

0

Percentage entropy margin for SS, ΔSS [%]

C

D
***

***

Groups in a quadrant [%]Homogeneous HOI parameters, θ k

10
15

*

10

θ5 0

5
0

100
80
60
40
20
0

−5
−10

−10
10

−15
3

4

5

6

θ3

The order, k

0

θ3+ - + - + - + -θ4+ + - - + + - -θ5+ + + + - - - -

−10
−10

0
10

θ4

Figure 4 | The groups that express HOIs exhibited alternating signs of homogenous HOIs at successive orders of interaction. (A) Scatter plots ofentropy margins for hHOIs versus entropy margins for SS. The same color indicates groups selected from the same slice culture. Filled circles indicate groups
H2 {Hss
exhibiting HOIs (DHOI . 3%). Dashed lines represent different proportions of homogenous HOI entropy margin explained by the SS term, b~.
H2 {HhHOI
We excluded 9 outliers from the plots. (B) Cumulative distribution function of this proportion, b, for the groups exhibiting HOIs (DHOI . 3%). (C) The
homogeneous HOI parameters up to the 6th order of the hHOI model. Each box covers 25th to 75th percentile, and whiskers represent 1.5 times the distancefrom the 25th to 75th percentile. Dots are outliers. The distributions at the 3rd, 4th, and 5th order deviated significantly from zero (two tailed sign test, *** and* represent significance level 0.001 and 0.05, respectively). (D) (Left) 3-dimensional plots of the homogeneous HOIs of the hHOI model. Outliers with
elements larger than 10 or less than 210 were excluded. (Right) A histogram of the number of groups that fell in 8 quadrants of the 
h3 ,
h4 ,
h5 parameter
space. The same color marks the same slice culture. The dotted horizontal line is the chance level (12.5%) with random HOIs.
SCIENTIFIC REPORTS | 5 : 9821 | DOI: 10.1038/srep09821

6www.nature.com/scientificreports

B

5
Cross correlogram

Neurons

Probability density

A

4
3

30
20
10

2

0
-1

1

-0.5

0

0.5

1

Time [s]

0

C

100

200

300
Time [s]

400

400

500

D

Total 450 samples

Number of groups

4%

300

SS groups
Other groups

200
100

0 0.2 0.4 0.6 0.8 1
Correlation coefficient

600

96%

0

Fraction of SS groups [%]0

20

Control
Picrotoxin

15
10
5
0

0

5

10

15

20

25

Percentage entropy margin for HOIs [%]30

200 400 600 800 1000
Bin size [ms]Figure 5 | Blocking inhibitory networks by PTX eliminated HOIs. The panels retain the same presentation format as in Fig. 1A, D and Fig. 3B. (A)Ensemble activity of 76 neurons in the CA3 area of a single hippocampal slice under bath application of PTX. (B) Distribution of correlation coefficients
between two event sequences (resolution of 100 ms) of all the pairs of neurons in 450 groups from 9 slices. Inset shows an average cross-correlogram. (C)Distribution of percentage entropy margin for HOIs. The groups that showed improved fitting with the SS term are marked in red (the SS groups, 4%).
Others (96%) are in gray. (D) The number of the SS groups with respect to the bin size for the analysis.

(DHOI . 3%). In these groups, the single SS term explained 80% of theentropy for homogeneous HOIs (interquartile range 48–92%). From
this result, we conclude that SS constitutes the dominant structure of
the homogeneous HOIs.
We next directly visualize the structure of homogeneous HOIs.
Figure 4C displays distributions of the homogeneous HOI parameters, hk of the hHOI models. (We only show the parameters for
k~3,4,5,6 although we fitted hHOIs up to the 10th order). The
homogeneous HOI parameters up to the fifth order but not higher
were significantly different from zero(two tailed
 sign test). The set ofthe homogeneous HOI parameters, h3 , 
h4 ,
h5 , from each group fellin a particular quadrant in the 3-dimensional space (negative triplewise, positive quadruple-wise, and negative quintuple-wise interactions, Fig. 4D), exhibiting an obviously biased direction. (If the set ofhomogeneous HOI parameters randomly fell in any quadrant, theprobability that the observed number of groups (, 54% of thegroups) would fall in any single quadrant would be less than10{15 ). Thus the structured homogeneous HOIs found up to the
5th order contributed to the excess SS found in 68% of the groups
with DHOI w3%. These results demonstrate that the structured HOIs
with alternating signs are an attribute of excess SS in local networksof hippocampal neurons.

inset) than that in the control (cf. Fig. 1C), much shorter than the
400 ms time window used to analyze the control condition.
Nonetheless, we used the same window-size, 400 ms, to test for thedeviation of SS from the pairwise model, except in Fig. 6D, where we
explored the dependency on bin sizes. Table 1 summarizes activity
rates, correlation coefficients, and probabilities of SS computed using
the 400 ms bin under control and PTX conditions. The average
probability of SS under the PTX conditions was much larger thanthat under the control condition. However, this frequent SS is
expected from the high pairwise correlation coefficients observed
under the PTX conditions. Indeed, the entropy explained by HOIs
was greatly diminished in the PTX data, indicating that the pairwise
model adequately explained population activity in almost all groups
under blockade of inhibition (Fig. 5C). Accordingly, the percentage
of groups that exhibited significant SS beyond the pairwise model
was considerably reduced from 16% down to 4% (Fig. 5C, red). The
considerable reduction of SS groups was observed whenever thewindow size larger than 200 ms was used in order to thoroughly
cover the synchronous events (Fig. 5D). We thus concluded thatan inhibitory network is necessary for neurons to produce both
frequent SS and weak pairwise correlations; the conjunction ofboth can only be explained by HOIs.

Simultaneous silence relies on network inhibition. Several
different biological mechanisms may underlie the observed
structure of HOIs. One such mechanism may be the inhibitory
networks in the hippocampal CA3 area. To test this hypothesis, we
examined population activity under bath application of GABAA
receptor antagonist picrotoxin (PTX) (Fig. 5A). When fast GABAA
mediated inhibitory networks were blocked by PTX, activities ofobserved neurons nearly completely synchronized with each other(Fig. 5B). The cross-correlogram exhibited a sharper peak (Fig. 5B

Simultaneous silence emerges in a population of thresholding
units that receive correlated input. Finally, we demonstrated thata simple model of neural population reproduces the structured HOIs
with alternating signs with respect to different orders of interaction
observed in the spontaneous activity of hippocampal neurons under
the control conditions. A population model known as the
Dichotomized Gaussian (DG) model33–36 simulates a population ofneurons that receive correlated Gaussian inputs, where each neuron
produces a binary output in response to its input by simple

SCIENTIFIC REPORTS | 5 : 9821 | DOI: 10.1038/srep09821

7www.nature.com/scientificreports
thresholding (Fig. 6A, See Methods). Despite the substantial
simplification, the spiking mechanism is similar to the one
assumed in networks of balanced excitatory and inhibitory
neurons: the mean input to each neuron is typically smaller than
the threshold and, therefore, spikes are induced by fluctuations in
the input. The DG model has been reported to reproduce neuralactivity patterns better than the pairwise model22. Figure 6B
displays simulated DG models using different strengths of input
correlations (Fig. 6B), including one that produces output
correlations similar to those found in experimental observations(see Table 1). The population exhibited asynchronous spiking
activity. We numerically computed the HOIs of the DG model for
N~10 (See Methods). The HOIs showed clear alternation in signs
with respect to the successive orders of interaction (Fig. 6C) anddemonstrated excess SS (Fig. 6D). These results show that theexperimentally observed SS with structured HOIs can arise fromthe conjunction of two ubiquitous biological features, i.e.,correlated input and spiking nonlinearity. Further, we demonstratethat SS can contain rich information of inputs provided to theobserved population of neurons. Figure 6E compares the signal-tonoise ratio for estimating the input correlation based on specificfeatures of population activities – activity rates and pairwiseAB‘an event’u1

cin=0.2, cout =0.05

cin=0.4, cout=0.13

100
80

0

cin=0.4

Neurons

correlated Gaussian
random variables−5u3

Discussion
We investigated the structure of HOIs in spontaneous activity ofneurons in the CA3 area of organotypic hippocampal slice cultures.
Most groups (, 70%) of neurons that expressed significant HOIs
(DHOI w3%) also exhibited excess SS (Fig. 3), and SS alone could
account for , 20% of the entropy explained by HOIs in these groups.
This result predicts significantly biased homogeneous HOIs withalternating signs at successive orders of interaction and our data
analysis confirmed this prediction (Fig. 4). We also found that SS
explained 80% of the entropy due to structured homogeneous HOIs.

0
−5

u2

correlations (Rate 1 Pair), SS in addition to activity rates andcorrelations (Rate 1 Pair 1 SS), and joint activity rates of allorders (Full). The signal-to-noise ratio, which is also called thelinear Fisher information37–39, quantifies the accuracy of estimatinga small change in an input parameter by an optimal liner decoder (seeMethods). The result shows that measuring SS in addition to activity
rates and correlations can provide nearly full information available
from the observation of all statistics. The results were qualitatively
the same when an input mean (or a threshold level) was estimatedinstead of input correlation.

0

60
40
20

−5
0

20

40

60

80

100

0

500

Time bins

0

500

1000

Time steps

E
1.4

1
0
−1

(cin=0)cout=0cout=0.05 (cin=0.2)
cout=0.13 (cin=0.4)

−2
−3

1.5

1.2

Signal-to-noise ratio

HOIs of the DG model

D

N=10

2

A parameter for SS, θ0

C

1000

Time steps

1.0
0.8
0.6
0.4

1.0
Full
0.5

Rate+Pair+SS
Rate+Pair

0.2
0

−4
1

2

3

4

5

6

7

8

9 10

The order of interactions, k

0.1

0.2

0.3

0.4

0.5

Input correlation coefficient, cin

0

0.1

0.2

0.3

0.4

0.5

Input correlation coefficient, cin

Figure 6 | The ensemble activity simulated by the Dichotomized Gaussian (DG) model exhibits alternating signs of HOIs depending on successiveorders of interaction. (A) Illustration of a DG model of 3 neurons. The traces in each panel represent correlated input variables, u~ðu1 ,u2 ,u3 Þ’ at different
time steps. The inputs are sampled from a multivariate Gaussian distribution, u*N ðc,LÞ, with mean vector c and covariance L (See Methods). Here weassume that the mean vector contains the same scalar element, {1:75, in order to yield the activity probability 0:04, a value close to the empiricallyobserved average activity rate (0.039 per 400 ms window). The off-diagonal elements of L are all fixed at cin ~0:4. The vertical lines above 0 in each panel
mark the time steps at which an input crosses the threshold. (B) Simulation of the DG model with 100 neurons with weak (Left) and strong (Right) input
correlations. The weak input correlation (cin ~0:2) in the Left panel yields a weak correlation coefficient (cout ~0:05) of output binary variables, whereas
the strong input correlation (cin ~0:4) in the Right panel yields a stronger correlation coefficient (cout ~0:13) of output binary variables. (C) The HOIs of asmall population N~10 from the DG model shows clear alternation of signs as the order of interaction k increases except for cin ~0. Negative interactions
occur at odd k and positive interactions occur at even k. (D) The parameter of SS, h0 , as a function of input correlation coefficient, cin . The dots marked incolor represent h0 at cin ~0, 0.2, and 0.4. (E) The signal-to-noise ratio of the input correlation coefficient, cin , as a function of cin in the population activityof the DG model (The solid purple line). The dotted lines are signal-to-ratios in the subset features of the population activity (blue: the activity rates ofindividual and pairwise neurons; orange: the activity rates of individual and pairwise neurons plus the SS rate).
SCIENTIFIC REPORTS | 5 : 9821 | DOI: 10.1038/srep09821

8www.nature.com/scientificreports
SS was robustly observed across a range of time-bins (Fig. 5D) andsizes of neural populations (Fig. 3C). Moreover, SS and the resulting
structure of HOIs arise in the simplest model of a neural populationthat possesses a spiking nonlinearity and correlated inputs (the DG
model), where additional observation of SS is sufficient to decode
most information of input conveyed by different orders of HOIs.
These results suggest that excess SS is an important and ubiquitous
characteristic of neural population activity that summarizes its lowdimensional structure in the combinatorial space of HOIs. We
identified alternating signs of HOIs up to the 5th order withstatistical significance in the analyzed data (Fig. 4C and D), andthe DG model displayed the predicted structure up to the highest
order of interaction (Fig. 6C). Based on these observations, we speculate that the predicted structure of HOIs beyond the 5th order should
be identifiable in future, given longer experimental recordings. We
also speculate that appropriate models5,15,21 of neural population
activity implicitly include excess SS as well as the resulting structureof HOIs, as demonstrated in the DG model.
Multiple biological mechanisms may underlie the high SS probability observed and alternating signs of HOIs. While we have shownthat even simple thresholding units with correlated inputs can reproduce this structure, we do not exclude contributions from othermechanisms. Indeed, we demonstrated the involvement of inhibitory
input in generating SS (Fig. 5). Under the blockade of GABAA receptors, activities of neurons were almost completely synchronized.
Therefore HOIs of the population activity were significantly diminished. It is expected that neurons are almost fully synchronized andfire regularly if inhibition is removed40. However, it may require
additional neuronal mechanisms with slow dynamics41,42 to robustly
account for the sparse synchronous activity observed in the current
data sets. While this result may simply indicate that inhibition isnecessary to place a network of neurons in a fluctuation-driven
regime43,44 for them to be sensitive to correlated input (c.f. Fig. 6),
it may alternatively suggest the existence of clustered inhibitory inputthat simultaneously shuts down a group of local neurons and produces excess SS. Inhibitory interneurons in the hippocampus havediverging connections to principal neurons45 and show powerfulcontrol over timing and rhythms of their spiking activity46,47. Suchinhibitory circuits are ideally suited to implement a winner-take-alllike competition among groups of neurons, which are common inmodels of hippocampal circuits aiming to reproduce place fields48,49.
Similar operations of hippocampal inhibitory circuits have also beensuggested for cellular assemblies50 and memory consolidation51.
Thus the excess SS in spontaneous activity reported here might berelated to functions requiring sparse information representation witha small fraction of active neurons. It is therefore interesting to see if
the same experimental manipulation of inhibition that influences, forexample, the sparse place field representation, also influences SS
during spontaneous activity.
Our study demonstrates that excess SS explains a large fraction of
the variability caused by complex HOIs in neural populations.
Although it was previously reported that HOIs fitted to several
representative activity patterns explain occurrence probabilities ofother general patterns18, this study did not normalize the model
probability distribution because of the computational complexity
associated with the normalization step. As a drawback, it was previously unknown how much of the variability associated with HOIs
was explained by a small number of representative activity patterns.
In contrast, the entropy maximization approach we have taken was
suitable to evaluate these quantities. More generally, virtually all
previous studies of HOIs14,19,20 attempted to fit multiple model parameters to the data rather than to extract the most prominent featurein the space of HOIs. Note that the hHOI model fitted to the data
reproduces an observed histogram of the number of active neurons
in time bins (i.e., a population spike-count histogram). Thus thehHOI model is equivalent to the K-pairwise model proposed in
SCIENTIFIC REPORTS | 5 : 9821 | DOI: 10.1038/srep09821

Tkačik et al.20 although the two models utilize different features ofactivity patterns to represent homogenous HOIs. We have found that
SS can parsimoniously summarize 80% of the specific structure of
hHOIs. Furthermore, successive orders of interaction have alternating signs. This resulting structure extends the negative triple-wiseinteractions previously found in local (v300mm) populations of 3
neurons17.
In sum, we demonstrate that representing HOIs using ‘‘silence’’
provides a much more concise description than the canonicalrepresentation based on ‘‘activity’’. We conclude that significant SS
is a ubiquitous feature in neural population activity that expresses
apparently diverse HOIs across different orders.

Methods
Recording method. Hippocampal slice cultures were prepared from postnatal day 7
Wistar/ST rats (SLC) (either male or female). Entorhino-hippocampal stumps were
cultivated on membrane filters using 50% minimal essential medium, 25% Hanks’balanced salt solution, 25% horse serum, and antibiotics in a humidified incubator at
37uC in 5% CO2 and were used for experiments on days 7 to 14 in vitro. On
experimental days, slices were washed with oxygenated artificial cerebrospinal fluid(aCSF) consisting of (mM) 127 NaCl, 26 NaHCO3, 3.3 KCl, 1.24 KH2PO4, 1.2 MgSO4,
1.2 CaCl2, and 10 glucose and bubbled with 95% O2 and 5% CO2. They were thentransferred to a 35-mm dish filled with 2 ml of dye solution and incubated for 40 min
in a humidified incubator at 37uC in 5% CO2 with 0.0005% Oregon Green 488
BAPTA-1 (OGB-1) AM (Invitrogen), 0.01% Pluronic F-127 (Invitrogen), and 0.005%
Cremophor EL (Sigma-Aldrich). They were recovered in aCSF for. 30 min andmounted in a recording chamber at 32uC and perfused with aCSF at a rate of 1.5–
2.0 ml/min for . 15 min. Hippocampal CA3 pyramidal cell layer was imaged at
10 Hz using a Nipkow-disk confocal microscopy (CSU-X1; Yokogawa Electric), a
cooled CCD camera (iXonEM1 DV897; Andor Technology), an upright microscopewith a water-immersion objective lens (16 3, 0.8 numerical aperture, Nikon).
Fluorophores were excited at 488 nm with a laser diode and visualized with a 507-nm
long-pass emission filter. The recording lengths varied from 600 sec to 3300 sec(600 sec (n 5 9); 1200 sec (n 5 4); 310, 610, 700, 900, 1100, 1800, 3300 s (n 5 1)).
Picrotoxin was bath-applied at a concentration of 50 mM to 9 slices (600 sec (n 5 7)
and 350 s (n 5 2)). After identification of cell types, the regions of interest (ROIs)
were carefully placed onto the cell bodies. The fluorescence change (DF/F) was
calculated as DF/F 5 (Ft–F0)/F0, where Ft is the fluorescence intensity at time t, and F0
is the baseline averaged for 50 s before and after time t. For neurons, event times werereconstructed from the onsets of Ca21 transients3,27,28. The signals were then inspectedby eye to remove erroneously detected noise. The data is available online (http://gaya.
jp/data). Under the same condition as described above, membrane potentials werewhole-cell recorded at I 5 0 from pyramidal cells (n 5 7) visually identified under
infrared differential interference contrast microscopy. Patch pipettes (326 MV) werefilled with a solution consisting of (in mM) 120 K-gluconate, 10 KCl, 10 HEPES, 10
creatine phosphate, 4 MgATP, 0.3 Na2GTP, and 0.2 EGTA. The signal was digitized at
10 kHz and filtered with a band of 1–2000 Hz. Liquid junction potentials were notcorrected. Experiments were performed with the approval of the animal experimentethics committee at the University of Tokyo (approval No. P24-6) and according to
the University of Tokyo guidelines for the care and use of laboratory animals. Allefforts were made to minimize the animals’ suffering and the number of animals used.
Selection of groups of neighboring neurons. From each slice, we selected 50 distinct
overlapping groups, each consisting of 10 nearest-neighbor neurons, based on thefollowing procedure. In each slice, we estimated the density of spatial distribution of
the cells in the recorded area of CA3 by an optimized 2-dimensional kernel density
estimation method52. We then sampled a spatial point from the estimated density,and selected the 10 neurons nearest to the point. We repeated this procedure until we
obtained 50 distinct groups (we discarded groups of neurons if the exactly same groupof 10 was previously selected). The neurons with low activity rates (less than 0.01 Hz)
were excluded from this analysis. In addition, we changed the number of neurons in agroup from 3 to 14 to investigate effect of the group size (Fig. 3C). We sampled up to
50 groups per slice following the same sampling procedure described above. Finally,
we sampled non-overlapping groups of N~3,    ,14 neurons. Note that we cansample only a small number of groups from each slice if groups are stochasticallysampled by the above-mentioned method. Thus, we took the following procedure toefficiently select non-overlapping groups from each slice. First, we fitted a 2dimensional Gaussian density function to the spatial distribution of cells in each slice.
We then determined the first principle component, and scored the positions ofneurons along the first principle axis. We selected neurons that are nearest neighbors
in terms of this score as a group. To determine the number of groups sampled fromeach slice, we computed the maximum number of non-overlapping groups that canbe sampled from each of all slices. We used the smallest number of groups among
them to sample an equal number of groups from each slice.
Model fitting and a test of simultaneous silence. First, we fit to binary population
8,9,31
activity data (see Fig.
1F in Results) the pairwise,hX
i maximum entropy modelX:::
hi x i z
hij xi xj , where xi is a binary variable of 0 or 1.
p2 (x1 , ,xN )* exp
i

ivj

9

www.nature.com/scientificreportsHere the parameters of the model, hi ,hij , were fitted by a maximum likelihood
principle, hMLE ~ arg max lðhÞ, where lðhÞ is the log likelihood of the data under thefeatures achieved by a general optimal nonlinear decoder as assessed by the Fisher
information.

hmodel. The nonlinear fitting was performed using a custom convex optimization
program in Matlab. We then fit to the same data a maximum entropy model thataugments the pairwise model with a single term to account for the observed
probability of SS in addition (the SS model, see Eq. 1 in Results). The increase inlikelihood over the pairwise model seen after adding a parameter for the SS is related
to the reduction in entropy by 2ðlss {l2 Þ~2T ðH2 {Hss Þ, where lss and l2 are the log
likelihood of the data under the SS model and the pairwise model respectively, and T
is the number of observed patterns (bins). Under the null hypothesis of no such SS
term, the variability in the model estimation due to finite samples make the differencein log likelihood following a x2 -distribution with one degree of freedom asx2 ð1Þ11. The p-value of the observed likelihood increase was computed using2ðlss {l2 Þ~
this null distribution. The p-values were further corrected by the BenjaminiHochberg-Yekutieli multiple comparison correction method that is applicable todependent tests, using Matlab code written by Groppe et al.53,54 This method controls
the proportion of tests that incorrectly declare significant SS (the false discovery rate).
A dichotomized gaussian (DG) model. The DG model is a threshold neuron modelwith Gaussian input signals33–36. The binary output of the i-th neuron (i~1,    ,N) is
given by Xi ~1 if ui w0 or Xi ~0 if ui ƒ0, where u~ðu1 ,u2 ,    ,uN Þ’ is drawn from amultivariate Gaussian distribution with mean c~ðc1 ,c2 ,    ,cN Þ’ and a covariance~ ðc,LÞ. Note that ’ describes matrix (or vector)
matrix L whose diagonal is 1 as uNtranspose. Here we consider a homogenous neuron pool: the mean is all fixed at
ci ~{h and the off-diagonal elements of L are all fixed at cin . The probability that
individual output neurons are in an active state is given by g1 ~Wð{hÞ, where W is
the one-dimensional cumulative distribution function (CDF) of a zero-mean, unit
variance Gaussian distribution. The probability of simultaneous activity of 2 neurons
is given by g2 ~W2 ð{h,cin Þ, where W2 is the 2-dimensional Gaussian CDF with zeromeans, unit variances, and a off-diagonal correlation coefficient cin . The correlation
coefficient between 2 output neurons is given by cout ~g2 {g21 g1 ð1{g1 Þ.
The probability distribution of population activity has a simple analyticalexpression
in this model. Note that the correlated inputs can be written aspﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃui ~ 1{cin vi z cin e{h, where vi is a unit variance white Gaussian noisevi *N ð0,1Þ specific to each neuron, and e*N ð0,1Þ is an input noise that is common
across all neurons. The conditional probability of a single neuron spiking given thecommon input e is given by33

ð
pﬃﬃﬃﬃﬃ
2
1 ?
1
1 h{e cine{u =2 du~f ðeÞ~ pﬃﬃﬃﬃﬃ
1{erf pﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
: ð3Þ

2
1{c
p
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
2p 
2pﬃﬃﬃﬃﬃinh{e cin = 1{cinThe probability that exactly m neurons are active and N{m neurons are inactive is
given
 
pDG ðmÞ~E Nm f ðeÞm ð1{f ðeÞÞN{m ,ð4Þ
where the expectation is performed
with respect to the common input noise, e. Notethat the binomial factor Nm sums all possible combinations of population activity
patterns with m active neurons. In order to obtain the probability mass function for
the finite population size N, we numerically computed the above equation. On theother hand, the same population-count probabilities are described by
hXm
i
 
F ðmÞhk {y :ð5ÞpðmÞ~ Nm exp
k~1 k
 
Here Fk ðmÞ~ mk (k~1,    ,N) is the k-th order feature of the hHOI model, whichcounts all combinations of choosing k neurons out of m active neurons, and y is anormalization factor. Thus, by solving linear equations, log pðmÞ~ log pDG ðmÞ for
m~0,1,    ,N, we obtain the parameters, hk (k~1,    ,N), and the normalization
factor y.
We quantify the signal-to-noise ratio for estimating the input correlation, cin , basedon the population activity of the homogenous DG model as follows. A small change incin is inferred from a vector of observation, Fobs ðmÞ~ðFs1 ðmÞ,    ,Fsr ðmÞÞ’ , whereindices sj [f1,    ,N g ðj ¼ 1,    ,r Þ specify a subset of r features that are taken into
account for the inference. The signal for detecting the input correlation is given by
LE ½Fobs ðmÞ and the noise of the observation is quantified by Cov½Fobs ðmÞ, where E½: is=Lcin
the expectation and Cov½: is the r|r covariance matrix calculated using pDG ðmÞ
defined above. Together, the signal-to-noise ratio is given by
LE½Fobs ðmÞ’LE½Fobs ðmÞ 38,39
Cov½Fobs ðmÞ{1. In the paper, we specifically consider
Lcin
Lcin
three types of observations as Fobs ðmÞ: the full observation ðF1 ðmÞ,    ,FN ðmÞÞ’ , theactivity rates of individual and pairwise neurons ðF1 ðmÞ,F2 ðmÞÞ’ , and these activity
rates plus the SS rate ðF1 ðmÞ,F2 ðmÞ,dm,0 Þ’ , where di,j represents Kronecker’s delta.
Notably, when all the features ðF1 ðmÞ,    ,FN ðmÞÞ’ are observed, the above signal-tonoise ratio becomes equivalent to the Fisher information55 of the input correlation,
L2i.e., E { 2 log pDG ðmÞ , and thus upper-bounds the accuracy of unbiased
Lcin
estimators of cin based on the population activity. See56,57 for information in subset

SCIENTIFIC REPORTS | 5 : 9821 | DOI: 10.1038/srep09821

1. Salinas, E. & Sejnowski, T. J. Correlated neuronal activity and the flow of neural
information. Nat Rev Neurosci 2, 539–550 (2001).
2. Reyes, A. D. Synchrony-dependent propagation of firing rate in iteratively
constructed networks in vitro. Nat Neurosci 6, 593–599 (2003).
3. Takahashi, N., Sasaki, T., Matsumoto, W., Matsuki, N. & Ikegaya, Y. Circuit
topology for synchronizing neurons in spontaneously active networks. Proc Natl
Acad Sci USA 107, 10244–102499 (2010).
4. Vidne, M. et al. Modeling the impact of common noise inputs on the network
activity of retinal ganglion cells. J Comput Neurosci 33, 97–121 (2012).
5. Ko€ster, U., Sohl-Dickstein, J., Gray, C. M. & Olshausen, B. A. Modeling higherorder correlations within cortical microcolumns. PLoS Comput Biol 10, e1003684(2014).
6. De la Rocha, J., Doiron, B., Shea-brown, E. & Reyes, A. Correlation between neuralspike trains increases with firing rate. Nature 448, 802–807 (2007).
7. Pitkow, X. & Meister, M. Decorrelation and efficient coding by retinal ganglion
cells. Nat Neurosci 15, 628–635 (2012).
8. Schneidman, E., Berry, M. J., Segev, R. & Bialek, W. Weak pairwise correlations
imply strongly correlated network states in a neural population. Nature 440,
1007–1012 (2006).
9. Shlens, J. et al. The structure of multi-neuron firing patterns in primate retina.
J Neurosci 26, 8254–8266 (2006).
10. Martignon, L. et al. Neural coding: higher-order temporal patterns in theneurostatistics of cell assemblies. Neural Comput 12, 2621–2653 (2000).
11. Nakahara, H. & Amari, S. Information-geometric measure for neural spikes.
Neural Comput 14, 2269–2316 (2002).
12. Tang, A. et al. A maximum entropy model applied to spatial and temporalcorrelations from cortical networks in vitro. J Neurosci 28, 505–518 (2008).
13. Yu, S., Huang, D., Singer, W. & Nikolic, D. A small world of neuronal synchrony.
Cereb Cortex 18, 2891–2901 (2008).
14. Montani, F. et al. The impact of high-order interactions on the rate ofsynchronous discharge and information transmission in somatosensory cortex.
Philos Trans. A Math Phys Eng Sci 367, 3297–3310 (2009).
15. Montani, F., Phoka, E., Portesi, M. & Schultz, S. R. Statistical modelling of higherorder correlations in pools of neural activity. Physica A 392, 3066–3086 (2013).
16. Ohiorhenuan, I. E. et al. Sparse coding and high-order correlations in fine-scalecortical networks. Nature 466, 617–621 (2010).
17. Ohiorhenuan, I. E. & Victor, J. D. Information-geometric measure of 3-neuron
firing patterns characterizes scale-dependence in cortical networks. J Comput
Neurosci 30, 125–141 (2011).
18. Ganmor, E., Segev, R. & Schneidman, E. Sparse low-order interaction network
underlies a highly correlated and learnable neural population code. Proc Natl
Acad Sci USA 108, 9679–9684 (2011).
19. Tkačik, G. et al. The simplest maximum entropy model for collective behavior in aneural network. J Stat Mech Theor Exp P03011 (2013).
20. Tkačik, G. et al. Searching for collective behavior in a large network of sensory
neurons. PLoS Comput Biol 10, e1003408 (2014).
21. Santos, G. S., Gireesh, E. D., Plenz, D. & Nakahara, H. Hierarchical interaction
structure of neural activities in cortical slice cultures. J Neurosci 30, 8720–8733(2010).
22. Yu, S. et al. Higher-order interactions characterized in cortical activity. J Neurosci
31, 17514–17526 (2011).
23. Shimazaki, H., Amari, S., Brown, E. N. & Grün, S. State-space analysis of timevarying higher-order spike correlation for multiple neural spike train data. PLoS
Comput Biol 8, e1002385 (2012).
24. Levy, W. B. & Baxter, R. A. Energy efficient neural codes. Neural Comput 8,
531–543 (1996).
25. Okamoto, K. et al. Ex vivo cultured neuronal networks emit in vivo-likespontaneous activity. J Physiol Sci 64, 421–431 (2014).
26. Ikegaya, Y., Aaron, G. & Cossart, R. Synfire chains and cortical songs: Temporalmodules of cortical activity. Science 304, 559–564 (2004).
27. Sasaki, T., Matsuki, N. & Ikegaya, Y. Metastability of active CA3 networks.
J Neurosci 27, 517–528 (2007).
28. Takahashi, N., Sasaki, T., Usami, A., Matsuki, N. & Ikegaya, Y. Watching neuronal
circuit dynamics through functional multineuron calcium imaging (fMCI).
Neurosci Res 58, 219–225 (2007).
29. Mizuseki, K. & Buzsáki, G. Preconfigured, skewed distribution of firing rates in the
hippocampus and entorhinal cortex. Cell Rep. 4, 1010–1021 (2013).
30. Buzsáki, G. & Mizuseki, K. The log-dynamic brain: how skewed distributions
affect network operations. Nat Rev Neurosci 15, 264–278 (2014).
31. Cocco, S., Leibler, S. & Monasson, R. Neuronal couplings between retinal ganglion
cells inferred by efficient inverse statistical physics methods. Proc Natl Acad Sci
USA 106, 14058–14062 (2009).
32. Okun, M. et al. Population rate dynamics and multineuron firing patterns insensory cortex. J Neurosci 32, 17108–17119 (2012).
33. Amari, S., Nakahara, H., Wu, S. & Sakai, Y. Synchronous firing and higher-orderinteractions in neuron pool. Neural Comput 15, 127–142 (2003).
34. Bethge, M. & Berens, P. Near-maximum entropy models for binary neuralrepresentations of natural images. Adv NIPS 20, 97–104 (2008).

10www.nature.com/scientificreports
35. Macke, J. H., Opper, M. & Bethge, M. Common input explains higher-ordercorrelations and entropy in a simple model of neural population activity. Phys Rev
Lett 106, 208102 (2011).
36. Montangie, L. & Montani, F. Quantifying higher-order correlations in a neuronalpool. Physica A 421, 388–400 (2015).
37. Series, P., Latham, P. E. & Pouget, A. Tuning curve sharpening for orientation
selectivity: coding efficiency and the impact of correlations. Nat Neurosci 7,1129–1135 (2004).
38. Beck, J., Bejjanki, V. R. & Pouget, A. Insights from a simple expression for linear
fisher information in a recurrently connected population of spiking neurons.
Neural Comput 23, 1484–1502 (2011).
39. Moreno-Bote, R. et al. Information-limiting correlations. Nat Neurosci 17,
1410–1417 (2014).
40. Brunel, N. Dynamics of sparsely connected networks of excitatory and inhibitory
spiking neurons. J Comput Neurosci 8, 183–208 (2000).
41. Mongillo, G., Barak, O. & Tsodyks, M. Synaptic theory of working memory.
Science 319, 1543–1546 (2008).
42. Curto, C., Sakata, S., Marguet, S., Itskov, V. & Harris, K. D. A simple model ofcortical dynamics explains variability and state dependence of sensory responses
in urethane-anesthetized auditory cortex. J Neurosci 29, 10600–10612 (2009).
43. Van Vreeswijk, C. & Sompolinsky, H. Chaos in neuronal networks with balanced
excitatory and inhibitory activity. Science 274, 1724–1726 (1996).
44. Shadlen, M. N. & Newsome, W. T. The variable discharge of cortical neurons:
Implications for connectivity, computation, and information coding. J Neurosci
18, 3870–3896 (1998).
45. Freundl, T. F. and Buzsaki, G. Interneurons of the hippocampus. Hippocampus 6,
347–470 (1996).
46. Royer, S. et al. Control of timing, rate and bursts of hippocampal place cells bydendritic and somatic inhibition. Nat Neurosci 15, 769–75 (2012).
47. Bartos, M., Vida, I. & Jonas, P. Synaptic mechanisms of synchronized gamma
oscillations in inhibitory interneuron networks. Nat Rev Neurosci 8, 45–56 (2007).
48. Rolls, E. T., Stringer, S. M. & Elliot, T. Entorhinal cortex grid cells can map tohippocampal place cells by competitive learning. Network 17, 447–465 (2006).
49. Cheng, S. & Frank, L. M. The structure of networks that produce thetransformation from grid cells to place cells. Neuroscience 197, 293–306 (2011).
50. Klausberger, T. & Somogyi, P. Neuronal diversity and temporal dynamics: theunity of hippocampal circuit operations. Science 321, 53–57 (2008).
51. Alvarez, P. & Squire, L. R. Memory consolidation and the medial temporal lobe: asimple network model. Proc Natl Acad Sci USA 91, 7041–7045 (1994).
52. Shimazaki, H. & Shinomoto, S. Kernel bandwidth optimization in spike rateestimation. J Comput Neurosci 29, 171–182 (2010).
53. Benjamini, Y. & Yekutieli, D. The control of the false discovery rate in multiple
testing under dependency. 29, 1165–1188 (2001).

SCIENTIFIC REPORTS | 5 : 9821 | DOI: 10.1038/srep09821

54. Groppe, D. M., Urbach, T. P. & Kutas, M. Mass univariate analysis of event-related
brain potentials/fields I: a critical tutorial review. Psychophysiology 48, 1711–1725(2011).
55. Cover, T. M. & Thomas, J. A. Elements of Information Theory. (John Wiley & Sons,
Inc., 1991).
56. Oizumi, M., Okada, M., Amari, S. Information loss associated with imperfectobservation and mismatched decoding. Front Comput Neurosci. 5, 9 (2011).
57. Koyama, S. On the relation between encoding and decoding of neuronal spikes.
Neural Comput. 1425, 1408–1425 (2012).

Acknowledgments
We thank Takuya Sasaki for his help in providing picrotoxin data, and Charles Yokoyama
and Alexandra V Terashima for their critical reading of our manuscript. This work was
supported by RIKEN Brain Science Institute (HS and TT), Kakenhi 22115003 (YI), andHealth and Labour Sciences Research Grants for Research on Regulatory Science ofPharmaceuticals and Medical Devices from MHLW (iNCENS: iPS-Non ClinicalExperiments for Nervous System) (YI).

Author contributions
H.S.: Conception and design; Analysis and interpretation of data; Drafting and revising thearticle. K.S.: Drafting and revising the article. T.I.: Acquisition of data; Revising the article.
Y.I.: Acquisition of data; Drafting and revising the article. T.T.: Conception and design;
Analysis and interpretation of data; Drafting and revising the article.

Additional information
Supplementary Information accompanies this paper at http://www.nature.com/
scientificreports
Competing financial interests: The authors declare no competing financial interests.
How to cite this article: Shimazaki, H., Sadeghi, K., Ishikawa, T., Ikegaya, Y. & Toyoizumi,
T. Simultaneous silence organizes structured higher-order interactions in neural
populations. Sci. Rep. 5, 9821; DOI:10.1038/srep09821 (2015).
This work is licensed under a Creative Commons Attribution 4.0 International
License. The images or other third party material in this article are included in thearticle’s Creative Commons license, unless indicated otherwise in the credit line; if
the material is not included under the Creative Commons license, users will need
to obtain permission from the license holder in order to reproduce the material. To
view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

11

www.nature.com/scientificreports

OPEN

Erratum: Simultaneous silence
organizes structured higher-orderinteractions in neural populations
Hideaki Shimazaki, Kolia Sadeghi, Tomoe Ishikawa, Yuji Ikegaya & Taro Toyoizumi
Scientific Reports 5:9821; doi: 10.1038/srep09821; published online 28 April 2015; updated on 09 October
2015
This Article contains typographical errors.
In the Results section under subheading ‘Simultaneous silence and HOIs of hippocampal neurons’“The parameters {θi,θij} were adjusted to fit these statistics.”
should read:“The parameters {θi,θij} were adjusted to fit these statistics.”
In the Results section under subheading ‘Simultaneous silence is a ubiquitous feature of HOIs’“By expanding the SS term into the standard HOI-coordinates, we obtainθ 0 ∏ (1 − x i ) = θ 0 − θ 0 ∑ x i + θ 0 ∑ x i x j − θ 0 ∑ x i x j x k + θ 0
i

i

i< j

i< j< k∑ x ix jx kx l.”i< j< k < lshould read:“By expanding the SS term into the standard HOI-coordinates, we obtainθ 0 ∏ (1 − x i ) = θ 0 − θ 0 ∑ x i + θ 0 ∑ x i x j − θ 0 ∑ x i x j x k + θ 0
i

i

i< j

i< j< k∑ x ix jx kx l − .”i< j< k < lIn the Result section under subheading ‘Alternating signs of HOIs predicted by SS’, equation (2)10phHOI (x1,  , x10) ∼ exp ∑θ ix i + ∑ θi1i 2x i x + ∑ θ k ∑ x i1x ik ,1ix2
i1< i 2k = 3 i1<  < i k ishould read:10

phHOI (x1,  , x10) ∼ exp ∑θ ix i + ∑ θ i1i2x i1x i2 + ∑ θ k ∑ x i1x ik ,i1< i 2k = 3 i1<  < i k iScientific Reports | 5:14893 | DOI: 10.1038/srep14893

1www.nature.com/scientificreports/
In the Methods section under subheading ‘A dichotomized gaussian (DG) model’,“The binary output of the i-th neuron (i = 1 ,  , N ) is given by X i = 1 if u i > 0 or X i = 0 if u i ≤ 0,
where u = (u1, u 2,  , u N )′ is drawn from a multivariate Gaussian distribution with mean∼γ = (γ1, γ 2,  , γ N )′ and a covariance matrix Λ whose diagonal is 1 as uN (γ, Λ).”should read:“The binary output of the i-th neuron (i = 1 ,  , N ) is given by X i = 1 if u i > 0 or X i = 0 if u i ≤ 0,
where u = (u1, u 2,  , u N )′ is drawn from a multivariate Gaussian distribution with meanγ = (γ1, γ 2,  , γ N )′ and a covariance matrix Λ whose diagonal is 1 as u ∼ N (γ, Λ).”In addition,“The correlation coefficient between 2 output neurons is given by c out = η 2 − η12/ η1(1 − η1).”should read:“The correlation coefficient between 2 output neurons is given by c out = (η 2 − η12 ) / [η1(1 − η1) ].”In addition,“The probability that exactly m neurons are active and N-m neurons are inactive is given”
should read:“The probability that exactly m neurons are active and N-m neurons are inactive is given by”
In addition,“The signal for detecting the input correlation is given by ∂E[F obs (m) ]/∂c inand”

should read:“The signal for detecting the input correlation is given by ∂E[F obs (m) ] and”
∂c inThis work is licensed under a Creative Commons Attribution 4.0 International License. The
images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under theCreative Commons license, users will need to obtain permission from the license holder to reproduce
the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Scientific Reports | 5:14893 | DOI: 10.1038/srep14893

2