PHYSICAL REVIEW E 93, 062416 (2016)

Clustering of neural code words revealed by a first-order phase transition
Haiping Huang* and Taro Toyoizumi†
RIKEN Brain Science Institute, Wako-shi, Saitama 351-0198, Japan(Received 8 March 2016; published 28 June 2016)
A network of neurons in the central nervous system collectively represents information by its spiking activity
states. Typically observed states, i.e., code words, occupy only a limited portion of the state space due toconstraints imposed by network interactions. Geometrical organization of code words in the state space, criticalfor neural information processing, is poorly understood due to its high dimensionality. Here, we explore theorganization of neural code words using retinal data by computing the entropy of code words as a function ofHamming distance from a particular reference codeword. Specifically, we report that the retinal code words inthe state space are divided into multiple distinct clusters separated by entropy-gaps, and that this structure isshared with well-known associative memory networks in a recallable phase. Our analysis also elucidates a special
nature of the all-silent state. The all-silent state is surrounded by the densest cluster of code words and locatedwithin a reachable distance from most code words. This code-word space structure quantitatively predicts typicaldeviation of a state-trajectory from its initial state. Altogether, our findings reveal a non-trivial heterogeneous
structure of the code-word space that shapes information representation in a biological network.
DOI: 10.1103/PhysRevE.93.062416
I. INTRODUCTION

Recent advances in multi-electrode recording techniques
allow simultaneous measurements of neural activity from alarge population of interacting neurons [1,2]. A populationof neurons encodes various information by its collective
spiking activity patterns, namely, neural code words [3]. These
code words are passed and interpreted by a downstream
circuit for further information processing. Characterizing theorganization of the code words is therefore critical for our
understanding of neural coding.
To characterize the distribution of code words, a maximum
entropy model [4] with pairwise interaction terms has been
fitted to neuroscience data [5,6]. This model that fits the firsttwo moments of activity statistics was reported to characterizereal data well in small groups of neurons. Importantly, thesestudies also suggest that code words are restricted due to neuralinteractions within a small subset of the state space, namely, thespace composed of all possible combinations of each neuron’sbinary activity. However, the geometrical organization of codewords is not well understood.
Interestingly, the code words of the well-known Hopfield
network [7] are also restricted within a small subset ofstate space due to strong constraints imposed by interactions
between neurons. The state space of the Hopfield networkis organized into multiple basins of attraction [8], withwhich a simple Glauber dynamics [9] can recall one ofmemorized patterns hinted by a distorted initial pattern. Thisis the so-called associative memory [7,8]. Although both theneuroscience model described above and the Hopfield network
belong to the pairwise maximum entropy model, it remains
largely unknown if their code-word spaces, composed of allcode words, share common features. Recent investigation ofretinal activity data revealed multiple local energy minima(LEM) in a fitted maximum entropy model [10]. However, it*†physhuang@gmail.com
taro.toyoizumi@brain.riken.jp

2470-0045/2016/93(6)/062416(11)does not provide how neural code words are geometricallyorganized because demonstration of the code-word spacestructure entails consideration of all possible states.
The high dimensionality of the state space prevents an
exhaustive search except in small networks, and standard
dimensionality-reduction techniques can easily abolish underlying structure by neglecting many relevant dimensions.
Hence, an efficient new technique is in need to visualize
the neural code-word space. One insight is that distance
between code words is an important factor that constrainsneural dynamics—previous experiments have shown thatstate transitions are mostly restricted to neighboring codewords and nearby code words are known to encode similar
information [11,12]. Based on this observation, we propose
the distance-constrained statistical mechanics analysis [13,14]
to concisely characterize the code-word space structure basedon the distance from a reference code word. In particular,
we present an advanced mean-field framework that computes
the entropy of code-words as a function of Hamming distancefrom any reference state. By applying this technique to both the
Hopfield network and retinal data, we explore their code-word
space structures, i.e., whether code words are divided into
multiple clusters.

II. RESULTS
A. Distance-constrained statistical mechanics analysis

We first introduce a statistical mechanics framework tocharacterize code word organization in the state space.
Let σi be binary activity of neuron i (i = 1, . . . ,N) andσ = (σ1 , . . . ,σN )T be a state vector, representing population
activity of N neurons. Here σi = 1 indicates that neuron iis active and σi = −1 indicates that neuron i is silent. The
symbol T represents the transpose operation.
According to the maximum entropy principle [4], theactivity state follows the Boltzmann distribution P (σ ) ∝exp(−βE(σ )), where β is the inverse temperature or neuralreliability (β = 1 unless otherwise indicated) and the energy

062416-1©2016 American Physical Society

HAIPING HUANG AND TARO TOYOIZUMIPHYSICAL REVIEW E 93, 062416 (2016)

(a)

(b)

FIG. 1. Code word organization of standard Hopfield model (P = 3) with low spiking rate constraint of stored patterns. (a) A first-orderphase transition in Hamming distance d when the coupling field x is tuned. The reference pattern is the one with four spikes [see (b)]. The
inset shows that the first order transition vanishes in the high temperature regime (β = 0.2). (b) entropy per neuron as a function of Hamming
distance from a reference stored pattern with different spike-counts. The pattern with zero spike-count is named all-silent (AS) state. The
curves correspond to the low-d branch of the hysteresis loop [see (a)]. The inset shows a trivial entropy landscape identical for all referencesin high-temperature regime.

E(σ ) = −hT σ − 12 σ T Jσ . h denotes a spiking bias vector andJ a functional coupling matrix. Geometrical organization ofcode words is studied by introducing a modified probability
distribution P (σ ) ∝ exp(−βE(σ ) + xσ T σ ∗ ), where coupling
field x is introduced to control the overlap σ T σ ∗ between state
σ and reference one σ ∗ . This perturbed probability measure
gives the free energy per neuron defined by
1log
exp(−βE(σ ) + xσ T σ ∗ )
βN
σ

1
=−
ddq exp (−Nβf (,q)),
βN

f ≡−

(1)

where
βf (,q) ≡ β − xq − s(,q)(2)is the energy- and overlap-dependent free energy that characterizes the probability of states having
 energy N  and overlap N q, and s(,q) ≡ (1/N) log σ δ( − E(σ )/N)δ(q −σ T σ ∗ /N)) denotes entropy (log-number of states) per neuron
with energy N and overlap N q. If the system-size N is large,
the integral in Eq. (1) is typically dominated by a combination(,q) that minimizes f (,q), i.e., f ≈ min,q f (,q).
We compute  and q that minimize f (,q) by applying theBethe approximation [15] (see the Appendix). By recursively
solving the mean field equation, we estimate a local (or global)
minimum of the free energy and (,q) corresponding to thisminimum. Notably, these values of  and q characterize theenergy and overlap of typically observed states (namely codewords), respectively. Meanwhile, the entropy of code wordss(,q) can be also computed according to Eq. (2). We define
Hamming distance (N − σ T σ ∗ )/2 that counts how manyneurons have distinct activity in state σ and reference state σ ∗ .
The typical value of the overlap q can be transformed to thetypical value of Hamming distance per neuron d = (1 − q)/2.
In the following sections, we omit the  dependency of theentropy and report it as a function of d, i.e., s(d).

B. Clustering of code words in the Hopfield model

Using the mean field method, we first investigate thestructure of code-word space in the Hopfield network [8,16].
In this model, the coupling between neuron i and j isμ μconstructed as Jij = 1/N Pμ=1 ξi ξj for a network of N =60 neurons, where P = 3 random binary patterns (indexedby μ = 1, . . . ,P ) are stored. In each pattern, stored activity
μξi of neuron i takes +1 with probability r = 0.0338 and −1
with probability 1 − r. r is chosen to fit the activity level ofretinal neurons we study in the next section. Note that, in the
Hopfield model, the neurons have zero spiking bias parameters(h = 0). According to the previous section, we compute thetypical distance d as we increase x from −3 to 3, and thendecrease it from 3 to −3 [Fig. 1(a)]. More precisely, after
the convergence of the mean field equations at some x, wechange x by a small amount and restart iteration from theprevious fixed point (see the Appendix). The reference stateσ ∗ is set to one of the stored patterns. Remarkably, we find afirst-order phase transition of d, characterized by the hysteresis
loop [Fig. 1(a)]. As we decrease x from high to low values, thetypical distance suddenly jumps at around x = −1.08 fromd ≈ 0.16 to d ≈ 0.97, implying a non-trivial structure of thecode-word space.
In order to more directly visualize the non-trivial structureof the code-word space, we plot the entropy of code words
computed at various distance d away from each stored pattern.
Only the entropy values corresponding to the low-d branchof the hysteresis loop are shown in the figure. As shownin Fig. 1(b), each stored pattern has a dense core of codewords around itself, which discontinuously falls off at somedistance. This indicates that code words are organized into
multiple clusters, separated by non-code-word states (i.e.,gaps). Among three stored patterns, the all-silent (AS) state
has the largest core due to the low spiking rate constraint of
stored patterns (small r).
This clustering results from the attractor structure [16] in the
retrieval phase of the model. Within the hysteresis loop, there
are two local minima of the free energy [Eq. (2)] competingwith each other. Low-d minimum corresponds to the nearby

062416-2

CLUSTERING OF NEURAL CODE WORDS REVEALED BY A . . .

Code-word space(a)

PHYSICAL REVIEW E 93, 062416 (2016)

Code-word space(b)

Clustering

No clustering

E( )

E( )

FIG. 2. Schematic illustration of two hypotheses on the organization of neural code words in the state space. Each binary neural code word(σ ) has an energy E(σ ). (a) When neural interactions are strong, neural code words can be organized into multiple clusters in the state space.
(b) When the neural population is sufficiently noisy, a trivial structure (a single cluster of neural code words) is observed.

code words of stored patterns (ξ ), while high-d minimum
corresponds to nearby codewords of corresponding reversed
patterns (−ξ ). Thus each stored pattern has distinct entropy
landscape surrounding it. The code-word space clustering isnecessary for successful memory retrieval in the Hopfield
network. In fact, in a high temperature regime (non-recallable
phase), the first-order transition and the non-trivial entropy
landscape are absent, as observed in the insets of Fig. 1. All
reference patterns display the same entropy landscape without
entropy gaps, and thus the patterns can not be distinguishedfrom each other. In this non-recallable phase of the Hopfield
network, there do exist multiple LEM [see Fig. 2(b)] undergreedy descent dynamics (GDD, see the Appendix) on theenergy surface, while the code-word space structure is trivial
without entropy gaps.

C. Clustering of code words in the retinal network

The next important question is how code words of a real
neural population are organized. To elucidate this question, we
analyze spiking activity data of populations of retinal ganglion
cells under a repeated naturalistic movie stimulus [10,17].
Although multiple LEM were previously found using thisdata set [10], it is still unknown if the observed network
has clustering of code words or not [Fig. 2(a)]. We therefore
characterize the geometrical organization of retinal code words
by applying the same method as used in the Hopfield network.
The neural spike trains in a population of N neurons
are binned with a 20 ms temporal resolution to have N dimensional spiking states σ . Spiking bias h and functional
coupling J are fitted to the spike train data to reproduce themean activity and pairwise correlation of the data (see the
Appendix). We choose randomly a network sample of the size
N = 60 from the neural data (the behavior reported below does
not change qualitatively when another sample is chosen, seeFig. 4). Despite no clear similarity in the connectivity structure
to the Hopfiled network, the retinal network displays thefirst-order phase transition with a hysteresis loop, qualitatively
resembling the Hopfield model [Fig. 3(a)]. This establishes
that code words of the retinal network are also clustered. Furthermore, by constructing an independent maximum entropy
model, where only the mean activity is fitted to the data with
J = 0, we show that the first-order phase transition disappears,

indicating that it is the non-trivial neural correlations that shape
the clustering of code words.
Figure 3(b) shows the entropy as a function of d whenneural code words of different spike-counts are selected as
references. Again, only the entropy values corresponding to
the low-d branch of the hysteresis loop are shown in the
figure. The high-d branch is not biologically plausible, since
the neural code is sparse. The entropy landscape is strongly
dependent of the reference. In general, the higher spike-counts
a neural code word has, the larger distance its entropy curve
extends over, enhancing the ability of the high spike-countcode word to come back to the sparse coding regime aroundthe AS state. To quantify this property, Fig. 3(c) plots the
maximum distance dmax at which the low-d branch in thehysteresis loop terminates as a function of spike-counts of
the reference code word. We find that dmax increases linearlywith the spike-counts (distance from the AS state) and theestimated slope is 0.9462 ± 0.0372. The slope close to one isalso observed in another typical example (see Fig. 4). Note thatthe distance to the AS state is typically smaller than dmax . Thisimplies that, even if the neural code word is far away from the
AS state, it still has easy access to the sparse coding regime
around the AS state within reasonable time, which highlights
the potential role of the AS state [18].
The AS state plays a special role here because the entropy
curve from the AS reference state grows much more rapidly
as a function of distance than from the other code words[Fig. 3(d)]. Indeed, its growth is close to the upper boundNgiven by the random code word limit [sub (d) = N1 ln(Nd)], inwhich every state is equally likely. This indicates that the AS
state has the densest core of code words around it, which
would facilitate frequent visits from other neural code words(see Fig. 5). As previously observed [10], a large portion ofneural patterns (about 94.25% of 2000 patterns) are observed
to evolve to the AS state by following GDD (see the Appendix).
Figure 3(d) reports the distance-dependent entropy landscape for some reference LEM code words (e.g., state a)obtained by running the GDD method starting from corresponding reference non-LEM code words (respectively statea  ) [see the corresponding multidimensional scaling (MDS)
map of LEM in Fig. 6(a)]. The result shows that each reference
has a different landscape, and at small d, the entropy arounda non-LEM code word is typically smaller than that for the062416-3

HAIPING HUANG AND TARO TOYOIZUMIPHYSICAL REVIEW E 93, 062416 (2016)

(a)

(b)

(c)

(d)

FIG. 3. Code word organization of the neural data (N = 60). (a) A first-order phase transition in Hamming distance when the coupling
field is tuned. The reference is a code word of ten spike-counts. The neural code-word space structure is shaped by the correlations in theneural spiking activity. The first-order transition disappears for independent model (IND). dmax defines the distance at which the low-d branchin the hysteresis loop terminates. (b) Distance-dependent entropy landscape from reference neural code words of different spike-counts. (c)Maximum distance dmax versus spike-counts of the reference states (distance from AS state). Five references for each spike-count are randomly
chosen. The line is a linear fit. (d) Distance entropy from neural code words (state a  ) and their corresponding LEM (state a). The corresponding
LEM are identified by GDD. The random code word limit is the upper bound.

corresponding LEM codeword. Moreover, for some states(e.g., 3 and 3 ), there exist two continuous parts separated
by a gap in the distance entropy curve. We shall elaborate
this phenomenon in the following section by studying a larger
population, where the effect becomes much more evident. Thisshows another clear evidence for the clustering of neural codewords.
D. More complicated structure observed
for large neural populations

The property of the neural codeword shown above is stillpreserved when large populations of neurons are considered.
In Fig. 7, we show the theoretical result computed on a
network sample of N = 100. As the network size grows, thenumber of LEM detected by GDD method also increases.
Accordingly, the internal structure of the code words becomes
more complicated [a rough visualization is given by the MDS
map, see Fig. 6(b)]. As shown in Fig. 7(a), there exist two
hysteresis loops separated by another monostable branch (twocurves for increasing and decreasing x coincide with eachother). These two successive hysteresis loops naturally arise if
there exist three deep minima in the energy landscape, where

sweeping x shifts a dominant contribution from one to another.
Fig. 7(b) shows that dmax grows with spike-counts (distancefrom the AS state). The growth is likely nonlinear in thiscase, perhaps induced by the complexity of the state space.
The fraction of neural code words that can reach the AS state
without in-between gaps reduces from the result of the previous
section to about 76.05%. Note that this number is still dominantcompared to the reachability of other detected LEM. Again,
the AS state has the densest surrounding core, characterized by
the rapid growth of the entropy with distance [Fig. 7(c)]. The
entropy landscape surrounding the AS state does not have asecond monostable branch beyond the first entropy gap, except
at a biologically implausible distance close to 1. This might bebecause there is no deep enough minima around the AS state.
In contrast, the entropy landscape surrounding some other reference code words, e.g., state 1 , exhibits a second monostablebranch beyond the first entropy gap [see Fig. 7(a)], likely
indicating that there is another deep minimum around them.
To demonstrate the implication of the entropy landscape,we study how distance from a local energy minimumchanges with time when the neural system explores thestate space. We use the local dynamics rule characterized by
the transition probability w(σi → −σi |Hi ) = e−2σi Hi , where

062416-4

CLUSTERING OF NEURAL CODE WORDS REVEALED BY A . . .

(a)

PHYSICAL REVIEW E 93, 062416 (2016)

(b)

(c)

(d)

(e)

FIG. 4. Entropy landscape of the neural data (N = 60, another typical example). (a) A first-order phase transition in Hamming distancewhen the coupling field is tuned. The transition disappears for an independent model (IND). (b) Distance entropy from reference neural codewords of different spike-counts. (c) Maximum distance dmax at which the low-d branch in the hysteresis loop terminates versus spike-counts of
the reference (distance from AS state). Five references for each spike-count are considered. The line is a linear fit (slope=1.041 ± 0.040). (d)Distance entropy from neural code words and their corresponding LEM. (e) Low-dimensional representation of LEM corresponding to (d) bymultidimensional scaling (MDS) analysis (y,x serve as coordinates). State 0 is the AS state.

Hi = hi + j Jij σj denotes the effective spiking bias of neuron i. Under this dynamics, states are sampled from the original
distribution P (σ ) ∝ exp (−βE(σ )). Note that the GDD rule
to obtain LEM allows only monotonically decreasing energy
on the energy surface. In contrast, the current dynamics rule
allows the energy to increase occasionally. Sampled distancefrom a reference local energy minimum σ ∗ is denoted byd0 (t) = (1 − σ T (t)σ ∗ /N)/2 where t denotes the time step.
The mean field prediction d0MF of a typical code-word distanceis given by setting x = 0 and initializing the iteration equation(see the Appendix) at σ ∗ . Note that x = 0 corresponds to thecase without distance-constraint, and thus takes into accountall code words in the cluster that σ ∗ belongs to. As expected,
this calculation predicts the fluctuation plateau of d0 (t) close
to the reference, as shown in Figs. 8(a) and 8(b). Note thatthe local dynamics escapes fast from the AS state [see theinset of Fig. 8(a)], which may be related to its very small core[Fig. 7(c)]. The same qualitative behavior holds for the smaller
network (N = 60, see Fig. 9) and when the neural dynamics
is simulated starting from a non-LEM code word.

062416-5

HAIPING HUANG AND TARO TOYOIZUMIPHYSICAL REVIEW E 93, 062416 (2016)

(a)

(b)

FIG. 5. Distance or spike-counts
shown in Fig. 3). The time-dependent Hamming
 evolution of the neural data (N = 60, typical example
N− i σi (t)σi (t+1)
,
and
the
time-dependent
spike-countsr(t)≡distance is defined as d(t) ≡i δ(σi (t) − 1). (a) The profile for four repeats.
2
Trial-to-trial variability is observed. (b) The profile for only one repeat. The AS state is frequently visited, and the neural network seems to
explore the state space by local moves.

(a)

(b)

5×10−5−5×10−5

FIG. 6. Low-dimensional representation of LEM (identified by GDD) by multidimensional scaling (MDS) analysis (y,x serve as
coordinates). State 0 is the AS state. (a) MDS map for a population of 60 neurons (see Fig. 3). (b) MDS map for a population of 100 neurons(see Fig. 7), in which 29 LEM are identified by GDD method. It becomes difficult to represent faithfully these LEM in a low-dimensional space(some information are lost), nevertheless, the map still shows how they are distributed.
062416-6

CLUSTERING OF NEURAL CODE WORDS REVEALED BY A . . .

(a)

(b)

(c)

FIG. 7. Code word organization of the neural data (N = 100).
(a) Two hysteresis loops are observed. The reference is state 1
[see (c)]. (b) Maximum distance dmax at which the low-d branchin the hysteresis loop terminates versus spike-counts of the reference(distance from AS state). Five references for each spike-count are
randomly chosen. (c) Entropy curve for different code words andtheir assigned LEM. Complex structure is observed for state 1, 1 ,and 2 .

III. DISCUSSION

In this work, we have established the resemblance ofcode word organization between the retinal network and the

PHYSICAL REVIEW E 93, 062416 (2016)

Hopfield network. In previous studies, the memory retrieval
function of Hopfield network was empirically compared tothe behavior of real networks [8,10,16,19]. However, notheoretical framework was proposed to build a solid relationship between these artificial and real biological networks.
In fact, they are naturally distinct in terms of detailed
parameters. Surprisingly, we have found that the two networks
both similarly organize their code words. The clustering ofcodewords has been identified by the first-order phase transition in the code-word distance. This transition is accompanied
by hysteresis loops, which becomes increasingly complex
as the network size grows. We have also revealed that the
AS state has a distinct role from other codewords. The
number of code words surrounding the AS state always grows
much more rapidly as a function of distance compared tothat surrounding other code words. Interestingly, despite thepresence of entropy gaps, most codewords even far away fromthe AS state could still have easy access to it because of theirsurrounding dense cores of code words typically extending
beyond the AS state. Thus, the most frequently observed AS
state plays a key role in serving as a hub facilitating neural
exploration of the code-word space.
The only knowledge a neuronal population can have comesfrom the population activity of interacting neurons. As shownin our study, there exists well-designed structure of code words
in the neural state space. The code words are partitioned into
multiple clusters separated by entropy gaps. Moreover, thisemergent property remains even if one-fourth of our data
is used to learn the model (see Fig. 10). Thus the revealedorganization structure is most likely an intrinsic property of
the retinal network, and downstream brain areas may benefit
from this structure for decoding purpose.
The clustering is functionally advantageous and intimately
related to the network function, i.e., pattern completion (errorcorrection) and pattern separation (discrimination ability).
Upon repeated presentations of the same visual stimulus,
the neural responses show strong trial-to-trial variability [10].
However, all code words belonging to the same core perhaps
encode the same feature of the semantic information [19]. This
property also allows the neural code to be robust against theubiquitous noise in nervous systems [20]. In an analogous
way to error-correcting codes [21,22], even if the neural code
word is corrupted by a small amount of noise, the dense core
structure still allows population coding of stimulus features.
Therefore, the non-trivial internal structure of the neuralcode-word space is useful for the neural population not only todiscriminate different neural activity patterns, but also to carryout error-correction [23,24].
The retina as an early visual system should adapt to thevisual stimulus distribution to efficiently transmit relevant
information to downstream brain areas. The energy landscape
shaped by the neural interactions likely depends on the naturalscene statistics [25]. It is therefore interesting to study theirrelationship under the current context.
The code-word space structure quantitatively predicts thefluctuation plateau of the simulated neural dynamics startingfrom LEM. Hence, our analytic framework establishes therelationship between the simulated neural dynamics andclustering of code words. In previous studies, the match
between spontaneous neural activity and the stimulus-evoked

062416-7

HAIPING HUANG AND TARO TOYOIZUMIPHYSICAL REVIEW E 93, 062416 (2016)

(a)

(b)

FIG. 8. Mean field theory predicts the plateau of the distance dynamics starting from LEM [d0 (t)]. (a) Typical trajectory observed forreference AS (inset), states 1 and 2 in simulations. The (solid, dashed, dotted) line is the theoretical prediction computed at x = 0 for eachreference. (b) The fluctuation plateau of d0 (t) is predicted by the mean field theory (d0MF ). Five trials from the same reference are considered
for each data point. Each trial lasts for 100 steps. Each step corresponds to N proposed flips. Note that in the inset of (a), one step corresponds
to one possible flip. LQ: lower quartile; MED: median; UQ: upper quartile.

activity increases during development especially for natural
stimuli [26], and the spontaneous activity outlines the regimeof evoked neural responses [27]. Our analysis might further

reveal how spontaneous neural activity is related to the
vocabulary of neural code words a neural circuit learns tointernally represent external worlds.

(a)

(b)

(c)

(d)

FIG. 9. Neural dynamics starting from LEM (d0 (t)). (a,b) the network with 60 neurons in the Fig. 3. (c,d) the network with 60 neurons
corresponding to Fig. 4. (a,c) typical trajectory observed in simulations. The (solid, dashed, dotted) line is the theoretical prediction computed
at x = 0 for each reference. (b,d) the fluctuation plateau of d0 (t) is predicted by the mean field theory (d0MF ). Five trials from the same reference
are considered for each data point. Each trial lasts for 100 steps. LQ: lower quartile; MED: median; UQ: upper quartile.
062416-8

CLUSTERING OF NEURAL CODE WORDS REVEALED BY A . . .

PHYSICAL REVIEW E 93, 062416 (2016)

ing one code word σ is written as P (σ ) ∝ exp (−E(σ )), where
hi σi −Jij σi σj .
(A1)E(σ ) = −iFIG. 10. Entropy landscape for N = 60. The same reference
codeword has ten spikes but two data of different lengths are learned.
The problem structure is not affected by the finite sampling of thedata.

Overall, our study provides an important step to understandthe stationary distribution of neural spiking patterns and its
functional relevance, which also sheds light on future studiesof the sensory processing in other brain areas.
ACKNOWLEDGMENTS

We are grateful to Michael J. Berry for sharing us theretinal data and letting us know recent similar results of hisgroup [25]. This work was supported by the program for Brain
Mapping by Integrated Neurotechnologies for Disease Studies(Brain/MINDS) from Japan Agency for Medical Research anddevelopment, AMED.
APPENDIX: METHODS
1. Simultaneous recordings of neural activity in populationsof retinal ganglion cells

The spiking activity of 160 retinal ganglion cells wascollected from a 450×450 μm patch of the salamander retina,
when a repeated naturalistic movie was presented. The visual
stimulus consists of 297 repeats of a 19s long movie clip beinga gray movie of swimming fish and swaying water plants ina tank (data courtesy of Michael J. Berry II, see experimental
details in the original paper [10,17]). The spike train data arebinned with the bin size τ = 20 ms reflecting the temporalcorrelation time scale, yielding about 280×103 binary neuralcodewords for model analysis.
2. Maximum entropy model

For a neuronal population of size N, the neural spike trainsof duration T are binned at temporal resolution τ , producing
M = T /τ samples of N -dimensional binary neural codewords. We use σi = +1 to indicate spiking activity of neuroni, and σi = −1 for silent activity. The neural responses torepeated stimulus are highly variable (so-called trial-to-trial
variability, see Fig. 5). To model the neural code-word statistics, we assign each code word σ a cost function (energy instatistical physics jargon) E(σ ), then the probability of observ-

i<j

The spiking bias hi and neuronal coupling Jij are constructedfrom the spike train data such that the spiking rate mi and the
pairwise correlation Cij under the model match those computed from the data. High energy state σ corresponds to low
probability of observation. This is a low-dimensional representation of the original high-dimensional neural code words,
since we need only N + N (N − 1)/2 model parameters.
To find the model parameters, we apply the maximum
likelihood learning principle corresponding to maximizing
the log-likelihood P (σ ) with respect to the parameters. The
learning equation is given by
ht+1= hti + η( σi data − σi model ),
i

(A2a)

Jijt+1 = Jijt + η( σi σj data − σi σj model ),(A2b)where t and η denote the learning step and learning rate,
respectively. The maximum likelihood learning shown here
has a simple interpretation of minimizing the Kullback-Leibler
divergence between the empirical probability and the model
probability [28,29]. In the learning equation [Eq. (A2)], thedata dependent terms can be easily computed from the binned
neural data. However, the model expectation of the spiking rate(magnetization in statistical physics) and pairwise correlationis quite hard to evaluate without any approximations. Here we
propose the mean field method to tackle this difficulty.
The statistical properties of the model [Eq. (A1)] can be analyzed by the cavity method in the mean field theory [30]. The
self-consistent equations are written in the form of message
passing (detailed derivation is given in Refs. [15,22]) as⎞⎛mi→a = tanh ⎝hi +(A3a)tanh−1 m̂b→i ⎠,b∈∂i\am̂b→i = tanh

mj →b ,b

(A3b)

j ∈∂b\iwhere ∂b\i denotes the member of interaction b except i, and∂i\a denotes the interaction set i is involved in with a removed.
a ≡ Jij and a ≡ (ij ). mi→a is interpreted as the message
passing from the neuron i to the interaction a it participatesin, while m̂b→i is interpreted as the message passing from theinteraction b to its member i. Iteration of the message passing
equation on the inferred model would converge to a fixed point
corresponding to a global (local) minimum of the free energy(in the cavity method approximation [15])
ln Zi +(|∂a| − 1) ln Za , (A4)F ≡ − ln Z = −iawhere Z is the normalization constant (partition function)of the model probability P (σ ). The free energy contribution of one neuron is − ln Zi = − ln x=±1 Hi (x), whereHi (x) ≡ exhi b∈∂i cosh b (1 + x m̂b→i ), and the free energy
contribution of one interaction is − ln Za = − ln cosh a −ln (1 + tanh a i∈∂a mi→a ). At the same time, the model

062416-9

HAIPING HUANG AND TARO TOYOIZUMIPHYSICAL REVIEW E 93, 062416 (2016)

spiking rate and multi-neuron correlation can also be estimatedas
mi = tanh hi +tanh−1 m̂b→i ,
(A5a)
b∈∂i

Ca =

tanh a +1 + tanh ai∈∂a mi→ai∈∂a mi→a

(A5b)

.

We have defined mi = σi and Ca =i∈∂a σi . Note thatthe iteration converges in a few steps at each learning stage,and estimated magnetizations as well as correlations are usedin the gradient ascent learning step. Here the multi-neuron
correlation is calculated directly from the cavity method
approximation [31] and expected to be accurate enough forcurrent neural data analysis. Another advantage is the low
computational cost. A more accurate expression could bederived from linear response theory [32] with much more
expensive computational cost.
Finally, one can also estimate the entropy of the model fromthe fixed point of the 
message passing equation. The entropy
is defined as S = − σ P (σ ) ln P (σ ), and it measures thecapacity of the neural population for information transmission.
More obvious variability of the neural responses implies larger
entropy value. Based on the standard thermodynamic relation,
S = −F + E, where E is the energy of the neural populationand given by
Ei +(|∂a| − 1) Ea ,
(A6a)
E=−
a

i

Ei =

hi



x=±1 xHi (x) +

x=±1 Gi (x)x=±1 Hi (x)

(A6b)

,

tanh a + i∈∂a mi→a,
1 + tanh a i∈∂a mi→a⎡exhi ⎣ b sinh b (1 + x m̂b→i )Gi (x) =
Ea =

(A6c)

a

b∈∂i

+x

⎤
b coshb (1 − tanh2

mj →b ⎦b)j ∈∂b\i×cosha (1 + x m̂a→i ).

(A6d)a∈∂i\b

3. Distance-constrained entropy analysis

To uncover the internal structure of the neural code-word
space, we introduce a modified probability measure [14]
⎛
⎞



1
βhi σi +βJij σi σj + xσi∗ σi ⎠,P (σ ) = exp ⎝
Z
i
i<j
i
(A7)
where β is the inverse temperature or neural reliability, and thecoupling field x is introduced to control the overlap between
the neural code word σ and a reference one σ ∗ .
The partition function Z can be approximated by a saddle
point analysis, i.e., Z exp(N s(d,) − βN  + xN q), fromwhich the free energy per neuron f (density) is given by−βf = s(d,) − β + xq, where  is the energy density(E/N), s(d,) the entropy density (S/N ) and q the typicalvalue of the overlap (σ T σ ∗ /N). Note that the Hamming
distance per neuron is related to the overlap by d = (1 − q)/2.
According to the double Legendre transform, the entropy
density is calculated via s(d,) = −βf + β − xq. eNs(d,)counts the number of valid configurations around the reference
satisfying both the distance constraint (d) and the energy
density (). Here β controls the energy level and x selects
the overlap
 or Hamming distance. The overlap q is given by
q = N1 i σi∗ mi with mi being calculated under the modified
probability measure. (x,β) obeys the following equations:∂s(d,)/∂d = 2x and ∂s(d,)/∂ = β. In this setting, theabove iteration equations [Eq. (A3)] remain unchanged exceptthat the bias is changed to hi → βhi + xσi∗ and the coupling is
rescaled as Jij → βJij . For the real neuronal network, the neural reliability β = 1, since the constructed biases and couplings
reflect the neural noise observed in the spike train data. For the
Hopfield model, higher β implies weaker thermal fluctuationand may correspond to a retrieval phase for pattern completion.
Note that to compute the entropy curve for metastable orunstable branches of distance-coupling field curve, one has tofix d by searching for compatible coupling field x, e.g., by the
secant method [33].
4. Finding local energy minima from neural activity pattern

To search for a local energy minimum starting from any
given neural activity pattern, we use greedy descent dynamics(GDD) in the energy landscape [10]. To be more precise, foreach neuron, we flip its activity if the flip will decrease theenergy. If we could not decrease the energy by flipping any
neuron’s activity, then a local energy minimum is identified.
Such minima are also called single-flip stable attractors,i.e., their energy cannot be decreased by flipping any single
neuron’s activity. We choose randomly a pattern set of size
2000 from the neural data to ensure that any two patterns are
rarely identical. By applying the GDD method, we identify a
LEM set whose size is much smaller than that of the patternset, with a large portion of patterns evolving to the all-silent
state. The number of LEM increases with the network size.
These LEM are then expressed in a low-dimensional space(called multidimensional scaling analysis (MDS) [34]).
MDS represents the proximity between LEM in the highdimensional space with some degree of fidelity by the distance
between points in the low dimensional space.
5. Independent maximum entropy model

In the case of fitting only the first moments (mean spiking
activity), the distance entropy can be computed exactly. The
result is given by1 
s(q(x)) =ln 2 cosh(hi + xσi∗ )
N i
−

1 
(hi + xσi∗ ) tanh(hi + xσi∗ ), (A8)N iiwhere hi = 12 ln 1+m
and q(x) = N1
1−mi062416-10



∗∗i σi tanh(hi + xσi ).

CLUSTERING OF NEURAL CODE WORDS REVEALED BY A . . .
[1] I. H. Stevenson and K. P. Kording, How advances in neuralrecording affect data analysis, Nat. Rev. Neurosci. 14, 139(2011).
[2] R. Yuste, From the neuron doctrine to neural networks, Nat. Rev.
Neurosci. 16, 487 (2015).
[3] F. Rieke, D. Warland, R. R. de Ruyter van Steveninck, and
W. Bialek, Spikes: Exploring the Neural Code (MIT Press,
Cambridge, 1997).
[4] E. T. Jaynes, Information theory and statistical mechanics,
Phys. Rev. 106, 620 (1957).
[5] E. Schneidman, M. J. Berry, R. Segev, and W. Bialek, Weakpairwise correlations imply strongly correlated network statesin a neural population, Nature 440, 1007 (2006).
[6] A. Tang, D. Jackson, J. Hobbs, W. Chen, J. L. Smith, H. Patel,
A. Prieto, D. Petrusca, M. I. Grivich, A. Sher, P. Hottowy, W.
Dabrowski, A. M. Litke, and J. M. Beggs, A maximum entropy
model applied to spatial and temporal correlations from cortical
networks in vitro, J. Neurosci. 28, 505 (2008).
[7] J. J. Hopfield, Neural networks and physical systems withemergent collective computational abilities, Proc. Natl. Acad.
Sci. USA 79, 2554 (1982).
[8] Daniel J. Amit, H. Gutfreund, and H. Sompolinsky, Spinglass models of neural networks, Phys. Rev. A 32, 1007
(1985).
[9] R. J. Glauber, Time dependent statistics of the ising model,
J. Math. Phys. 4, 294 (1963).
[10] G. Tkacik, O. Marre, D. Amodei, E. Schneidman, W. Bialek,
and M. J. Berry, Searching for collective behavior in a large
network of sensory neurons, PLoS Comput. Biol. 10, e1003408(2014).
[11] M. Stopfer, V. Jayaraman, and G. Laurent, Intensity versus identity coding in an olfactory system, Neuron 39, 991(2003).
[12] M. M. Churchland et al., Stimulus onset quenches neuralvariability: A widespread cortical phenomenon, Nat. Neurosci.
13, 369 (2010).
[13] L. Dall’Asta, A. Ramezanpour, and R. Zecchina, Entropy
landscape and non-Gibbs solutions in constraint satisfaction
problems, Phys. Rev. E 77, 031118 (2008).
[14] H. Huang, K. Y. M. Wong, and Y. Kabashima, Entropy landscapeof solutions in the binary perceptron problem, J. Phys. A: Math.
Theor. 46, 375002 (2013).
[15] M. Mézard and A. Montanari, Information, Physics, andComputation (Oxford University Press, Oxford, 2009).
[16] D. J. Amit, H. Gutfreund, and H. Sompolinsky, Information
storage in neural networks with low levels of activity, Phys.
Rev. A 35, 2293 (1987).

PHYSICAL REVIEW E 93, 062416 (2016)[17] O. Marre, D. Amodei, N. Deshmukh, K. Sadeghi, F. Soo, T. E.
Holy, and M. J. Berry, Mapping a complete neural populationin the retina, J. Neurosci. 32, 14859 (2012).
[18] H. Shimazaki, K. Sadeghi, T. Ishikawa, Y. Ikegaya, and T.
Toyoizumi, Simultaneous silence organizes structured higherorder interactions in neural populations, Sci. Rep. 5, 9821
(2015).
[19] E. Ganmor, R. Segev, and E. Schneidman, A thesaurus for aneural population code, eLife 4, e06134 (2015).
[20] A. A. Faisal, L. P. J. Selen, and D. M. Wolpert, Noise in the
nervous system, Nat. Rev. Neurosci. 9, 292 (2008).
[21] N. Sourlas, Spin-glass models as error-correcting codes, Nature
339, 693 (1989).
[22] H. Huang and H. Zhou, Cavity approach to the Sourlas code
system, Phys. Rev. E 80, 056113 (2009).
[23] T. Hosoya, S. A. Baccus, and M. Meister, Dynamic predictive
coding by the retina, Nature 436, 71 (2005).
[24] S. Sreenivasan and I. Fiete, Grid cells generate an analog errorcorrecting code for singularly precise neural computation, Nat.
Neurosci. 14, 1330 (2011).
[25] A. Loback, J. Prentice, M. Ioffe, and M. Berry, Investigating thestructure of the retinal ganglion cell population response probability landscape, Computational and Systems Neuroscience(cosyne), 2016.
[26] P. Berkes, G. Orban, M. Lengyel, and J. Fiser, Spontaneous
cortical activity reveals hallmarks of an optimal internal modelof the environment, Science 331, 83 (2011).
[27] A. Luczak, B. Peter, and K. D. Harris, Spontaneous events
outline the realm of possible sensory responses in neocortical
populations, Neuron 62, 413 (2009).
[28] S. Cocco and R. Monasson, Adaptive cluster expansion for the
inverse Ising problem: Convergence, algorithm and tests, J. Stat.
Phys. 147, 252 (2012).
[29] H. Huang, Sparse Hopfield network reconstruction with 1
regularization, Eur. Phys. J. B 86, 484 (2013).
[30] M. Mézard and G. Parisi, The Bethe lattice spin glass revisited,
Eur. Phys. J. B 20, 217 (2001).
[31] A. Braunstein, A. Ramezanpour, R. Zecchina, and P. Zhang,
Inference and learning in sparse systems with multiple states,
Phys. Rev. E 83, 056114 (2011).
[32] H. Huang and H. Zhou, Counting solutions from finite samplings, Phys. Rev. E 85, 026118 (2012).
[33] J. Nocedal and S. Wright, Numerical Optimization (Springer,
Berlin, 2006).
[34] M. A. A. Cox and T. F. Cox, Multidimensional scaling, in
Handbook of Data Visualization, Springer Handbooks Comp.
Statistics (Springer, Berlin/Heidelberg, 2008), pp. 315–347.

062416-11