Recurrent infomax generates cell assemblies, neuronalavalanches, and simple cell-like selectivity
Takuma Tanaka
ttakuma@mbs.med.kyoto-u.ac.jp
Department of Morphological Brain Science,
Graduate School of Medicine, Kyoto University, Japan.
Takeshi Kaneko
Department of Morphological Brain Science,
Graduate School of Medicine, Kyoto University, Japan.
CREST, JST.
Toshio Aoyagi
Department of Applied Analysis and Complex Dynamical Systems,
Graduate School of Informatics, Kyoto University, Japan.
CREST, JST.
July 28, 2008

Abstract

Recently, multineuronal recording has allowed us to observe patterned firings, syn1

chronization, oscillation, and global state transitions in the recurrent networks of central nervous systems. We propose a learning algorithm based on the process of information maximization in a recurrent network, which we call “recurrent infomax” (RI).
RI maximizes information retention and thereby minimizes information loss throughtime in a network. We find that feeding in external inputs consisting of information
obtained from photographs of natural scenes into an RI-based model of a recurrent
network results in the appearance of Gabor-like selectivity quite similar to that existing in simple cells of the primary visual cortex. We find that without external input,
this network exhibits cell assembly-like and synfire chain-like spontaneous activity as
well as a critical neuronal avalanche. In addition, we find that RI embeds externallyinput temporal firing patterns to the network so that it spontaneously reproduces thesepatterns after learning.

RI provides a simple framework to explain a wide range ofphenomena observed in in vivo and in vitro neuronal networks, and it will provide anovel understanding of experimental results for multineuronal activity and plasticity
from an information-theoretic point of view.

1

Introduction

Recent advances in multineuronal recording have allowed us to observe phenomena in the networksof the central nervous system (CNS) that are much more complex than previously thought to exist.
The existence of interesting types of neuronal activity, such as patterned firings, synchronization,
oscillation, and global state transitions has been revealed by multielectrode recording and calcium
imaging (Nadasdy et al., 1999; Cossart et al., 2003; Ikegaya et al., 2004; Fujisawa et al., 2006;
Sakurai and Takahashi, 2006). However, in contrast to the rapidly accumulating experimental2

data, theoretical works attempting to account for this wide range of data have been slower tomaterialize. These new data are partly explained by the classical hypotheses proposed purely on
theoretical grounds, such as the “cell assembly” of Donald Hebb (Hebb, 1949). However, to explain
a wider range of data, we have to extend the classical hypotheses on the basis of mathematics andinformation sciences.
We hypothesize that these characteristic types of neuronal activity in CNSs can be explainedby a theoretical model based on “infomax.” The process of information maximization (infomax(Linsker, 1988)) maximizes the information transmission from the input to the output of a feedforward network (Fig. 1A). In this paper, information is defined as in the information theory
proposed by Shannon 60 years ago (Shannon, 1948), which has been successfully applied to electrical engineering, computer science, and physics (Cover and Thomas, 2006). In information theory,“information” is measured on the basis of the probability P (x) that a system takes state x. For
example, P ([1, 0, 0, 1]) = 0.01 means that the relative frequency of occurrence that the first andfourth neurons fire and the second and third ones remain silent is one percent over the duration ofa long trial. Mutual information I(X; Y ) of two discrete random variables X and Y with a joint
probability distribution P (x, y), x ∈ X , y ∈ Y and marginal probability distribution P (x) andP (y) is defined by
I(X; Y ) =XX

P (x, y) logx∈X y∈Y

P (x, y)= I(Y ; X),
P (x)P (y)where X and Y are the sets of states. Taking logarithms to base 2, we can measure the mutual
information in bits. Mutual information I(X; Y ) is the information shared by input X and output
Y . In other words, it measures the reduction in the uncertainty of X due to the the knowledge
of Y and vice versa. Thus, maximizing the mutual information of the input and output improves
the information transmission in a feedforward network. It has been proposed that infomax infeedforward networks may provide an explanation of the stimulus selectivity of neurons in CNSs
3

(Tsukada et al., 1975; Atick, 1992; Bell and Sejnowski, 1995; Olshausen and Field, 1996; Bell andSejnowski, 1997; Lewicki, 2002). However, CNSs contain not only feedforward but also recurrent
synaptic connections (Fig. 1B), which endow networks with many interesting phenomena, someof which have been reported recently and several researchers have attempted to model (Diesmann
et al., 1999; Maass et al., 2002; Buonomano, 2005; Vogels and Abbott, 2005; Teramae and Fukai,2007). Therefore, we attempted to extend infomax to the case of recurrent networks, in which the
input to the neurons at time t consists of their own output at time t − 1 (Fig. 1C).
More specifically, a learning algorithm based on infomax in feedforward networks generates
information-efficient representation of the input in the output neurons of the feedforward network(Fig. 2A1 and A2). This algorithm adjusts the connection weights to realize the most efficient
information transfer from the input to the output. In this way, a network with small mutual information of input and output, that is, large information loss (Fig. 2A1), evolves into a network that
preserves a larger percentage of information (Fig. 2A2) through this algorithm. If the optimization
based on infomax is applied to a recurrent network in which the input to the neurons at time t
consists only of their own output at time t − 1, the mutual information of two successive states,b is maximized; that is, the information loss through time is minimized. We call this form
I(X; X),of infomax “recurrent infomax” (RI). An algorithm based on RI readjusts the connection weightsof the recurrent network to change a random network with large information loss (Fig. 2B1) into an
information-efficient network (Fig. 2B2). The role of RI is to allow a recurrent network to optimize
the synaptic connection weight in order to maximize information retention and thereby minimize
information loss by maximizing the mutual information of the temporally successive states of the
network.
In this paper, proposing a learning algorithm based on RI, we find that feeding in external
inputs consisting of information obtained from photographs of natural scenes into an RI-based

4

model of a recurrent network results in the appearance of Gabor-like selectivity quite similar tothat existing in simple cells of the primary visual cortex (V1). More importantly, we find thatwithout external input, this network exhibits cell assembly-like and synfire chain-like stereotyped
spontaneous activity (Hebb, 1949; Abeles, 1991; Diesmann et al., 1999) and a critical neuronalavalanche (Beggs and Plenz, 2003; Teramae and Fukai, 2007; Abbott and Rohrkemper, 2007). RI
provides a simple framework to explain a wide range of phenomena observed in in vivo and invitro neuronal networks, and it should provide a novel understanding of experimental results for
multineuronal activity and plasticity from an information-theoretic point of view.

2

Methods

Here we briefly describe our recurrent network model, leaving the details of the derivation toAppendix A. In this model, N neurons are connected according to the weight matrix Wij , andtheir firing states [xi (t) = 1 (fire) and 0 (quiescent)] at time step t are synchronously updated totime step t + 1. The firing state xi (t + 1) of neuron i at time step t + 1 is determined stochasticallywith the firing probabilitypi (t + 1) =

pmax
 P
,
1 + exp − j Wij (xj (t) − p̄j ) + hi (t)(1)

where hi (t) is the threshold of neuron i, and pmax is the maximal firing probability. When the
maximal firing probability pmax = 0.5, a neuron fires on average once every two time steps, even if
the neuron receives a sufficiently strong excitatory input at every time step. A small value of pmax
thus makes the firing of the neurons quite unreliable. In contrast, if pmax is close to 1, it is highly
probable that a strong input makes a neuron fire. Thus, pmax determines the reliability with whicha model neuron fires in response to an input.
To fix the mean firing probability of neurons i to p̄i , we update the threshold of neuron i, hi (t),5

at each step according tohi (t + 1) = hi (t) + ∆hi (t + 1)
= hi (t) + (xi (t + 1) − p̄i ),(2)

where the learning rate  for the threshold is set to 0.01 in all simulations. Eq. 2 fixes the mean
firing probability of neuron i in a manner that the threshold rises when the neuron fires and thethreshold falls when the neuron remains silent. When the firing states and the thresholds are
updated by Eqs. 1 and 2 for a sufficiently long sequence of time steps, hi (t) stops increasing or
decreasing and starts fluctuating around a certain value. Then, the time-average of the second
term of the righthand side of Eq. 2 vanishes, and thereby the time-average of xi (t), that is, thefiring rate of neuron i becomes equal to p̄i . Thus, the mean firing probability is fixed to p̄i .
Input xi (0) to the neurons at the first step t = 1 of the simulation was set to 0, and in thefollowing steps xi (t) was determined stochastically with Eq. 1. Unless otherwise stated, the neurons
in the model network do not have other inputs than their outputs at the previous step, and thereby
the dynamics of the network are completely determined by Eqs. 1 and 2 (Fig. 3A).
We performed simulations in blocks consisting of 20,000-100,000 time steps, updated Wij at
the end of each block, and then started the calculation for the next block (Fig. 3B). Outputs of
the neurons at the last step of block b − 1 were given as inputs to the neurons at the first step ofblock b. A simulation consists of 500-15,000 blocks.
To maximize information retention, our recurrent network starts from a random weight Wijinitialand develops toward an optimized network with Wijoptimized . The evolution of the weight matrix isdetermined by the gradient ascent algorithm,
Wij (b + 1) = Wij (b) + ∆Wij (b)= Wij (b) + η
6

∂I(b)
,
∂Wij

(3)

where Wij (b) is the connection weight Wij in block b and I(b) is the mutual information of two
successive states of the network in block b and η is the learning rate. To avoid Wij increasing
without bound, it is bounded above and below by wlimit and −wlimit , respectively.
We define the approximate mutual information in block b of two states separated by n − 1 stepsbyI (n) (b) = log |C| −1
log |D(n) |,
2





where
 E11 · · · E1N 


 .

.
.

..
.. 
C =  ..
,
EN 1 · · · EN N



(n)

D(n)

Eij

(n)

 E11 · · · E1N E1b1 · · · E1Nb 


 .

.
.
.
..
...

 ..
.
.
.
.
.
.
.






(n)
(n) 
E
 N 1 · · · EN N EN b1 · · · EN Nb 
,
= 


E (n) · · · E (n) E (n) · · · E (n) 
 b11
b
b
b 
b
1N
1b
1
1N


 .

.
.
.
.
.
 .
..
..
..
..
.. 
 .





(n)
(n)
(n)
(n)
ENb 1 · · · ENb N ENb b1 · · · ENb Nb
1 X=(xi (t) − p̄i )(xj (t) − p̄j ),Tt∈T (b)

(n)

Ebıj

=

1 X
(xi (t + n) − p̄i )(xj (t) − p̄j ),Tt∈T (b)

(n)

Eib

=

1 X
(xi (t) − p̄i )(xj (t + n) − p̄j ),Tt∈T (b)

(n)

Ebıb

1 X
=
(xi (t + n) − p̄i )(xj (t + n) − p̄j ),Tt∈T (b)T (b) is the set of latter half of steps in block b, and T is half of the number of steps containedin a block, that is, T = #T (b). The connection weights Wij are updated using correlation in thelatter half of steps in a block to let hi (t) converge in earlier half of steps in this block after Wij
7b to
was updated. I (1) is an approximation of mutual information of two successive states, I(X; X),be maximized (see Appendix A for derivation).
The gradient of the mutual information with respect to connection weight is approximated by
∂I(b)
≈
∂Wkl

1 X
(1 − δij )(Ebıbk Ebl + Ebıl Ebbk )(2(C −1 )ji − (D−1 )ji − (D−1 )j+N i+N )2 1≤i≤N1≤j≤N

12− ((1 − 2p̄k )(1 − 2p̄l )Ebkl + p̄k p̄l (1 − p̄k )(1 − p̄l ) − Ebkl)((D−1 )l k+N + (D−1 )k+N l ),
2
(1)

(1)

(1)

where Ebıj = Ebıj , Eib = Eib , Ebıb = Ebıb , and D = D(1) . Fig. 4 shows that the approximate mutual
information increased through this algorithm in the learning process of the simulation shown in
Fig. 5.
All models in this paper can be fully characterized by parameters N (50-432), p̄i (0.002-0.05),
pmax (0.5-0.95), η (0.2-20),  (0.01), and wlimit (100-1000). Parameter values used in simulations
are included in figure captions. At the beginning of the simulation, Wij was drawn from a uniform
distribution on [−0.5, 0.5] and hi was set to 0.
In the simulation depicted in Fig. 5, the input image was converted to a gray-scaled image andthen high-pass filtered using an exponential filter with the frequency responseL(f ) = exp(−(f /fc )4 ),where the cut-off frequency is fc = 256 (Olshausen and Field, 1997). The processing in the earlyvisual systems such as the retina and lateral geniculate nucleus can be regarded as high-passfiltering, and the output neurons correspond to the neurons in V1. The filtered image data wasused to generate firing patterns of input neurons by taking 12×12 randomly selected image patches
and then converting these to 288 binary inputs. The on-input and off-input neurons fired only when
the intensities of the corresponding pixels had positive and negative signs, respectively. For eachpixel i of the input, α|di | was compared to a random value u drawn from a uniform distribution on8[0, 1], where di is the intensity of pixel i = 1 . . . 144, and α is a constant parameter. If α|di | > u,
the state of the corresponding input neuron was set to 1, and if α|dn | ≤ u, it was set to 0. We
set the parameter α to fix the mean firing probability of the input neurons around 0.15. Under
this condition, the pixels caused the firing of the input neurons with a probability proportional toits intensity, except for 5% of the pixels, whose intensities were larger than 1/α. The simulation
program was written in C++.

3

Results

We first observed the behavior of this model network under external input. Image patches froma photograph preprocessed by a high-pass filter were used as the external input (Fig. 5A). The
neurons in this network were divided into three groups: 144 on-input and 144 off-input neurons,
and the 144 output neurons were randomly selected from the network (Fig. 5B1). Pixels withpositive and negative values in a randomly selected 12 × 12 image patch excited the corresponding
on-input and off-input neurons, respectively. The states of the input neurons were stochastically
set to 1 or 0 with firing probabilities proportional to the intensities of the corresponding pixels,
whereas the states of the output neurons were not set by the external input (see Methods fordetails). Instead, the firings of these neurons were determined by Eq. 1 with pmax = 0.95. Initially,
the connection weight Wij was a random matrix (Fig. 5C1), and we found that output neurons did
not exhibit clear selectivity with respect to the external input from the input neurons (Fig. 5D1)
upon averaging the image patches that evoked firings in an output neuron. After learning, however,the network self-organized a feedforward structure from the on-input and off-input neurons to the
output neurons (Fig. 5B2,C2). The output neuron became highly selective to Gabor functionlike stimuli (Fig. 5D2), exhibiting behavior quite similar to the selectivity of simple cells in the9V1 cortex (Hubel and Wiesel, 1959). Our optimization algorithm based on RI hence caused themodel network to become organized into a feedforward network containing simple cell-like output
neurons. It has been proven that the infomax accounts for the selectivity of simple cells (Belland Sejnowski, 1995, 1997). Bell and Sejnowski (1997) argued that the natural image patches arecomposed of independent localized edges such as Gabor functions and that these components canbe recovered by maximizing the mutual information of the input and the output. We thus see that
this result is consistent with the previous studies based on information theory.
In the simulation described above, the external input was fed into a network with high response
reliability (pmax = 0.95). Next, we examined the evolution of the spontaneous activity in a neuronal
network without external input. In this network, the approximate mutual information I (1) of
two successive states was maximized, and the approximate mutual information I (n) of two statesinterleaved with n − 1 steps after learning became larger than I (n) before learning (Fig. 6A). We
supposed that this improvement in information retention was a result of the emergence of repeatedactivity in the network. To identify repeated activity in the model network, we defined a repeatedpattern as a spatial pattern of neuronal firings that occurs at least twice in the latter half of atest block (Fig. 6B). Coloring repeated patterns consisting of ≥ 3 firing neurons in raster plotsof the network (Fig. 6D1,D2), we found that the number of repeated patterns increased after
learning. Several patterns were repeated in a sample of 250 steps as seen in Fig. 6D2, where therepeated patterns are indicated by consistently colored circles and connected by lines. Moreover,some patterns appeared to constitute repeated sequences. For example, sequence A, composedof the magenta, orange, and purple patterns, appears three times in Fig. 6D2. To quantify theincrease in repetition, we tabulated the numbers of occurrences of repeated patterns and sequences,and compared these numbers before and after learning (Fig. 6C). We found that both repeatedpatterns and repeated sequences increased significantly after learning. This indicates that the10present algorithm embeds not only repeated patterns but also repeated sequences of firings into
the network structure as a result of the optimization.
When a pattern in a sequence is activated at one step, it is highly probable that the next
pattern in that sequence will be activated at the next step. This predictability means that thestate of the network at one time step shares much information with the state at the next time step.
In contrast, when the dynamics of a network is highly stochastic and thereby repeated patterns
are rare, we cannot predict which pattern follows a given pattern nor reduce the uncertainty of
the next pattern by using the knowledge of the present pattern. In this case, mutual informationof two successive states is low. Sequences must be repeatedly activated and the network must be
deterministic in order to efficiently retain information in a recurrent network. Hence, we conclude
that the repeated activation of an embedded sequence is an efficient way to maximize information
retention in a recurrent network. These repeated patterns and sequences have been experimentallyobserved in vivo (Skaggs and McNaughton, 1996; Sakurai and Takahashi, 2006; Yao et al., 2007)
and in vitro (Cossart et al., 2003; Ikegaya et al., 2004), and their existence is suggested by thetheory of cell assemblies proposed by Hebb (1949) and the theory of synfire chains proposed byAbeles (1991). We thus see that RI accounts for the appearance of cell assemblies, sequences, andsynfire chains in neuronal networks.
In the simulations shown above, a small fraction of connections grew especially strong in the
network after learning (Fig. 6E2). So we ask the question, is the existence of a small number ofstrong connections a sufficient condition for the efficient information transfer? To answer this,
we randomly shuffled the components of the weight matrix of the network after learning shownin Fig. 6, and we found that shuffled networks exhibited lower mutual information and a smaller
number of occurrences of repeated sequences (Fig. 7, A and B). Thus, the existence of strong
connections does not necessarily imply that the network is efficient in retaining information. RI

11

improves information retention in recurrent networks while randomly introducing strong connections does not.
We next examined the behavior of the same spontaneous model in the case that the maximal firing probability was small (pmax = 0.5). For small pmax , the number of identically repeatedsequences is small, and the network seems to lose structured activity. However, we found characteristic network activity consisting of firing in bursts (Fig. 8A2), which are defined as consecutive
firing steps that are immediately preceded and followed by “silent” steps, with no firing. We foundthat after learning, the distribution P (s) of the burst size s, which is the total number of firings
in a burst, obeys a power-law distribution P (s) ∝ sγ with γ ≈ −1.5, whereas, before learning, wehave P (s) ∝ exp(−αs) (Fig. 8C). This result is consistent with experimental results. Recently,
Beggs and Plenz (2003) recorded the spontaneous activity of an organotypic culture from the cortex using multielectrode arrays. Defining an avalanche similarly to our bursts following a period
of inactivity, they found that the size distribution of avalanches is accurately fit by a power-law
distribution with exponent −1.5. To explain this, they argued that a neuronal network is tuned tominimize the information loss and that this is realized when one firing induces an average of one
firing at the next step. They showed that this condition yields the universal exponent −3/2, using
the self-organized criticality of the sandpile model (Bak et al., 1987; Harris, 1989). This condition
also holds for the present network, because, after learning, each neuron with pmax = 0.5 had twostrong input connections and two strong output connections on average (Fig. 8B2). The universal
exponent −3/2 was observed in the network for small pmax (Fig. 8C), but not for pmax = 0.95. Actually, the size distribution of bursts P (s) in the system did not exhibit a power-law distribution,and displayed several peaks, reflecting the existence of stereotyped sequences (data not shown).
We thus conclude that RI embeds information-efficient structures in which one firing induces on
average one firing at the next step in a network with small pmax .

12

To reveal the essential mechanism responsible for the behavior described above, we returned to
the recurrent network with an external input (Fig. 9). It has been observed that the hippocampal
firing sequences in the awake state are repeated during sleep (Skaggs and McNaughton, 1996;
Louie and Wilson, 2001) and that the spontaneous spiking activity in the visual cortex mimics
the movie-evoked response after repeated exposure to a movie (Yao et al., 2007). We investigatedwhether or not the firings presented during the learning period are replayed by the present model
after the learning. In the learning blocks, we repeatedly stimulated neurons 1, 3, and 2 in sequence(Fig. 9A1,B1). In the learning blocks, the state of neuron 1 was set to 1 (fire) at random intervals
ranging from 50 to 99 steps (time step t). At t + 2, the state of neuron 3 was set to 1, and at
t + 6, the state of neuron 2 was set to 1. In the successive test block, in which only neuron 1
was stimulated externally (Fig. 9A2), the firing of neuron 1 was followed by spontaneous firingsof neurons 3 and 2 (Fig. 9B2, arrows). In addition, the spontaneous firing of neuron 1 triggers
the sequence containing the firings of neurons 3 and 2 (Fig. 9B2, double arrows). The formof the weight matrix after learning reveals that a feedforward structure starting from neuron 1(1 → 7, 34 → 3, 5 → 49 → 18 → 11, 28 → 2) was embedded in the network (Fig. 9C). This
structure self-organizes in the network because, as we saw above, embedding a sequence of firings
into the network structure is an efficient way to retain information. It is thus seen that RI embeds
externally input temporal firing patterns into the network by producing feedforward structures,
and, as a result, the network can spontaneously reproduce the patterns.

4

Discussion

In this study, we have found that infomax in recurrent networks acts to optimize the network
structure by maximizing the information retained in the recurrent network. Many previous papers

13

concerning infomax in feedforward networks (Linsker, 1988; Atick, 1992; Bell and Sejnowski, 1995,
1997; Lewicki, 2002) have suggested that the stimulus selectivity of neurons in CNSs is accounted
for by infomax in feedforward networks. In contrast, although infomax in recurrent networks has
been studied, infomax is applied only to small recurrent networks that can be studied by using arandom search (Ay, 2002). This is because the analysis of recurrent networks is complicated byhistory-dependent dynamics due to the recurrent connections. In the present model, approximating
the mutual information of two successive states with second order correlations of neuronal firings,
we succeeded in deriving an algorithm that maximizes information retention in recurrent networks.
The present model reproduced the self-organization of simple cell-like selectivity shown in theprevious models and we successfully extended these previous results to the spontaneous activity
characteristic to recurrent networks. In the context of a simple maze task, for example, these
repeated patterns can be regarded as memory traces representing spatial cues and relationshipbetween successive items, and they have been supposed to help an animal in solving the mazetask (Dragoi and Buzsáki, 2006). An internal representation of the external input is essential inadaptation to environments, and the internal representation is constructed by RI in the form of
feedforward structures.
We have found that infomax in recurrent networks reproduces self-organization of cell assemblies
and neuronal avalanches. In contrast, most previous theoretical studies on cell assemblies, synfire
chains, and neuronal avalanches investigated the dynamics of neuronal firings on a network inwhich a feedforward structure underlying this characteristic type of activity had been embedded(Diesmann et al., 1999; Beggs and Plenz, 2003; Teramae and Fukai, 2007). Although these models
successfully reproduced experimental results, they could not explain how the embedded network
structure emerges. A recent theoretical study suggested that neuronal avalanches are accounted
for by a simple model for the growth of dendritic and axonal processes (Abbott and Rohrkemper,

14

2007). It seems that this model self-organizes a network structure which maximizes retained
information as in our model.
In our model, the network structure self-organized by the optimization algorithm resulted insimple cell-like activity, repeated sequences, and neuronal avalanches. Through evolution, animals
have acquired CNSs, which are extremely efficient information processing devices that improve an
animal’s adaptability to various environments. It is thus quite natural that these phenomena canbe regarded as a result of the optimization of information retention. Thus, in this paper and ourmodel, we have focused on information retention in a recurrent network although CNSs should be
optimized not only for information retention but also for categorization and generalization. On
the other hand, previous studies showed that synaptic plasticity rules experimentally observed andtheoretically proposed optimize the information transmission of individual synapses (Toyoizumi
et al., 2005; Pfister et al., 2006). Thus, neuronal networks with local plasticity rules optimized
to retain information could reproduce the experimental results of repeated activity patterns andavalanches. However, the learning rule of the present model is not local and requires global
information. We can optimize the activity of, for example, the half of the neurons in the network
if we approximate the mutual information of these N/2 neurons using the N/2 × N/2 correlation
matrix and update the connection weights among these neurons, leaving other connection weights
unchanged. Then, we observe that the occurrences of repeated sequences increases after thislearning but not as much as in the simulation shown in Fig. 6 (data not shown). Even though
this learning rule requires the information on only the half of the neurons in the network, thisrule is not local and requires global information on the activity of these N/2 neurons in thesystem.

To overcome this problem, our next goal is to derive a biologically plausible plasticity

rule in a bottom-up way employing RI, and to compare this rule with experimentally obtained
plasticity rules. We believe that RI will help us to understand the meaning of in vivo and in15

vitro experimental results, particularly to characterize the spontaneous activity of neurons in thecontext of information theory.

A

Algorithm

Here we describe the algorithm to maximize the mutual information of the present state, X, andb of the network.
the next state, X,
b = [xi (t+1)]
N neurons receive as input an output x = [xi (t)] at time t and generate an output x
at time t + 1. Neuron i takes two states, a firing state, xi = 1, and a non-firing state, xi = 0. The
firing probability of neuron i at time t + 1 is given by Eq. 1. We assume that Wij can take positive
and negative values, with positive and negative Wij corresponding to excitatory and inhibitory
connections, respectively. The threshold hi (t) evolves according to Eq. 2 and fixes the mean firing
probability of neuron i to p̄i .
To derive the algorithm that maximizes the mutual information of consecutive states, we first
approximate the entropy of the state, H(X), and the entropy of the joint distribution of twob Let P (x) be the probability that the state of the network is x = [xi ],successive states, H(X, X).
b) be the probability that the states of the network at consecutive steps are x andand P (x, xb = [bxxi ], respectively. Then, these entropies are defined by
H(X) = −X

P (x) log P (x),

x

b = −
H(X, X)Xb) log P (x, xb).
P (x, xx,bx

If the distribution of the state x is given by a Gaussian distribution with the correlation matrix E11 · · · E1N 


 .

.
.
..
.. 
C=
 ..
,




EN 1 · · · EN N
16

where Eij = E[(xi − p̄i )(xj − p̄j )], the entropy of the state is
H(X) =

N1
log |C| + (1 + log 2π)22(Cover and Thomas, 2006), and the entropy of the joint distribution of two successive states x andb is given byxb =H(X, X)1
log |D| + N (1 + log 2π)
2if this joint distribution is Gaussian with correlation matrix

 E11

 .
 ..



E
 N1
D=

E
 b11

 .
 .
 .


ENb 1

···

E1N

...

..
.

E1b1 · · ·
..
.

···

EN N EN b1

···

Eb1N

Eb1b1

...

..
.

..
.

···



ENb N ENb b1

E1Nb 

.. 
..
.
. 



· · · EN Nb ,· · · Eb1Nb 



.
..
.. 
.



· · · ENb Nb

where Ebıj = E[(bxi − p̄i )(xj − p̄j )], Eib = E[(xi − p̄i )(bxj − p̄j )], and Ebıb = E[(bxi − p̄i )(bxj − p̄j )].
b is given by
Therefore, the mutual information of two successive states x and x
b = H(X) + H(X)b − H(X, X)
b = log |C| − 1 log |D|.
I(X; X)
2(4)

b is identical to theHere we have used Ebıb = Eij , that is, the fact that correlation matrix of x
correlation matrix C of x. We assume that recurrent infomax is realized by maximizing the valueof the function in Eq. 4.
b are not Gaussian because of the discreteness of the
Although the distributions of x and x
neuronal states, this approximation gives a good estimate of the mutual information. We comparedb of two consecutive steps with this approximation. Fig. 10 shows
the mutual information I(X; X)
that mutual information is fit quite well by the form log |C|− 21 log |D|. Because this approximation17

b of Nrequires only correlation matrices, it enables us to estimate the mutual information I(X; X)
neurons, whose calculation in its original form requires the joint probability distribution of 22N
realizations of the firing states.
In addition, the quantity in Eq. 4 is a good index of the information retained in a recurrent
network even when it deviates significantly from the value of the mutual information. Maximizing
Eq. 4 results in the decorrelation of the state x due to log |C|, as well as in the increase ofb, owing to − 21 log |D|. A strong
the correlation between the the state x and the next state x
correlation between the states of the network at two successive steps increases the amount of
information transmitted over time, and strong decorrelation among the neurons at a step increases
the information capacity of the network. Thus, Eq. 4 is an effective measure of the information
retained in the recurrent network. Another advantage of using Eq. 4 as the value function is that
this function can be calculated by using only the second-order correlations. Although higher-ordercorrelations are useful in estimating the mutual information, calculating higher-order correlations istime consuming in numerical simulations and complicates the theoretical analysis. In the following
derivation of the algorithm, we use Eq. 4, and thus employ an approximation of the mutual
information in which the contribution of the higher-order correlations is not taken into account.
Similarly, we assume that there exists no correlation higher than second order among the xi
and x
bj . Then, we can assume that the probability of a state x is given by P (x) = z1 (x)/Z1 , where
!
Z1 =

X
x

z1 (x) =

X

exp

X

x

Jij (xi − p̄i )(xj − p̄j )ij

is the partition function, in which x runs over all realizations of the firing states. (Each summation
for which no range is expressed is assumed to run from 1 to N .) The variable Jij is dependent on
and determined by the second-order correlation matrix C according toX
1 X(xi − p̄i )(xj − p̄j ) exp
Jkl (xk − p̄k )(xl − p̄l )
Z1 x
kl
18

!
=

X
(xi − p̄i )(xj − p̄j )P (x) = Cij , (5)xand thereby it is dependent on Wkl . We assume that J is a symmetric matrix, that is, Jij = Jji ,without losing generality. Jij does not represent a real connection strength between neurons i andj, but rather the firing correlation between them. A positive and a negative Jij imply that thefirings of neurons i and j are positively and negatively correlated, respectively. In other words, we
assume that the state x is generated by a Boltzmann machine with connection strength Jij , andthat this Boltzmann machine has been trained to produce the correlation C (Hinton and Sejnowski,1983). We do not have to solve Eq. 5 to obtain the value of Jij , and as we see in the following,calculating the derivative of Jij with respect to Wkl suffices to maximize the value of the function
in Eq. 4. Next, we assume that conditional probability, P (bx|x), is given by z2 (bx|x)/Z2 (x), where!
Z2 (x) =Xz2 (bx|x) =

X
b
x

b
x

exp

X

Wij (bxi − p̄i )(xj − p̄j ) −

ij

X

hi (b
xi − p̄i ) .

(6)iAlthough we assume this for the general case, it is exactly correct in the case pmax = 1, because inthis case, from Eq. 6, we obtain
!
X
z2 (b
x|x)
1 Y
1 Y
=
exp
z3 (b
xi |x),Wij (bxi − p̄i )(xj − p̄j ) − hi (bxi − p̄i ) =P (bx|x) =Z2 (x)Z2 (x) iZ2 (x) ij
and we recover the firing probability of neuron i at time t + 1,
P (bxi = 1|x) =z3 (bxi = 1|x)= pi (t + 1),z3 (bxi = 0|x) + z3 (bxi = 1|x)where we have set the state of neuron i at step t to xi (t) = xi in Eq. 1. Hence, formulating thepartition function Z of the system in the form
Z =

X

=

X

b)
z(x, xx,bx

z1 (x)z2 (b
x|x)

x,b
x

!
=

X
x,b
x

exp

X
ij

Jij (xi − p̄i )(xj − p̄j ) +X

Wij (bxi − p̄i )(xj − p̄j ) −ijb) = z(x, xb)/Z.
we can write the joint probability P (x, x
19

X
i

hi (b
xi − p̄i ) ,Differentiating the correlation Eij with respect to Jkl , we obtain
∂EijZ
∂Jkl

∂E[(xi − p̄i )(xj − p̄j )]
∂Jkl
X
∂
b)
=
(xi − p̄i )(xj − p̄j )P (x, x
∂Jkl
x,b
x


∂ 1 Xb)
=(xi − p̄i )(xj − p̄j )z(x, x∂Jkl Z=x,bx

b)1 ∂Z X
1X∂z(x, xb)
− 2=
(xi − p̄i )(xj − p̄j )(xi − p̄i )(xj − p̄j )z(x, xZ∂Jkl
Z ∂Jkl
x,b
x

x,b
x

= E[(xi − p̄i )(xj − p̄j )(xk − p̄k )(xl − p̄l )] − E[(xi − p̄i )(xj − p̄j )]E[(xk − p̄k )(xl − p̄l )]= Eijkl − Eij Ekl ,where the superscript Z indicates that Eij is regarded as a function of the independent variables
Jkl , Wkl , and hk , although Jij and hi are dependent on Wkl . Preceding in the same way, we findthe following relations:∂EbıZb= Ebıbbkl − Ebıb Ebkl ,∂Wkl∂EbıZb
= −Ebıbbk + Ebıb Ebk = −Ebıbbk ,∂hk∂EiZ= Eikl − Ei Ekl = Eikl ,∂Jkl∂EbıZb= Ebıbkl − Ebıb Ekl ,∂Jkl∂EbıZj= Ebıj bkl − Ebıj Ebkl ,∂Wkl
where Ebıbbkl = E[(bxi − p̄i )(bxj − p̄j )(xk − p̄k )(xl − p̄l )], Ebıbbk = E[(bxi − p̄i )(bxj − p̄j )(bxk − p̄k )], Eikl =E[(xi − p̄i )(xk − p̄k )(xl − p̄l )], Ebıbkl = E[(bxi − p̄i )(bxj − p̄j )(xk − p̄k )(xl − p̄l )], and Ebıj bkl = E[(bxi −p̄i )(xj − p̄j )(bxk − p̄k )(xl − p̄l )], and we have used Ei = E[xi − p̄i ] = 0 and Ebı = E[bxi − p̄i ] = 0.
To obtain∂C, we have to calculate∂Wkl
X ∂EijZ ∂Jmn
∂Eij
=
,
∂Wkl
∂Jmn ∂Wkl
mn
20

because C and Eij are determined by Jmn according to Eq. 5, and, Jij is dependent on Wkl . Statesb obey the same distribution, and thereby the dependency of Jij on Wkl is determined byx and x∆Ebıb = ∆Eij ,∆Ei = 0.
Thus,

∂Jmn
is given by the solution of the following:∂Wkl∆Ebıb = ∆Eij =⇒X ∂EbıZbkl∂Jkl∆Jkl +X ∂EbıZb∂Wkl
kl∆Ei = 0 =⇒∆Wkl +k

X ∂E Zi∂Jklkl

X ∂EbıZb
∂hk

∆hk =

X ∂EijZkl∂Jkl∆Jkl ,∆Jkl = 0.

Rearranging terms, we obtain
X(Ebıbkl − Ebıb Ekl − Eijkl + Eij Ekl )∆Jkl −X

Ebıbbk ∆hk = −

k

kl

X

Eikl ∆Jkl = 0.

X
(Ebıbbkl − Ebıb Ebkl )∆Wkl ,

(7)

kl

(8)

kl

Because we have assumed that there exist no higher-order correlations, we substitute the higher-

21

order correlations in these equations with the second-order correlations as follows:p̄i (1 − p̄i )(1 − 2p̄i ) i = j = k,Eijl ≈0otherwisep̄i (1 − p̄i )(1 − 2p̄i ) bı = b
=bkEbıbbl ≈
,



0
otherwise(1 − 2p̄i )(1 − 2p̄k )Eiki = j, k = l





Eijkl − Eij Ekl ≈ (1 − 2p̄i )(1 − 2p̄j )Eij + p̄i p̄j (1 − p̄i )(1 − p̄j ) − E 2 i = k, j = l or i = l, j = k ,
ij








0
otherwise

Ebıj bkl − Ebıj Ebkl ≈Ebıbkl − Ebıb Ekl ≈(1 − 2p̄i )(1 − 2p̄j )Ebıj + p̄i p̄j (1 − p̄i )(1 − p̄j ) − Ebı2j bı = b
k, j = l0otherwise(1 − 2p̄i )(1 − 2p̄k )Ebık + p̄i p̄k (1 − p̄i )(1 − p̄k ) − Ebı2k bı = b
, k = l



0

,

,

otherwise

Ebıbbkl − Ebıb Ebkl ≈(1 − 2p̄i )(1 − 2p̄l )Ebıl bı = b
=b
k



E b Ebl + Ebıl E b
bık
b
k

.

otherwise

Here we have assumed that terms containing correlations among three or more variables are small
and can be set to zero, except in the last approximation. In the last approximation, we haveb isassumed Ebıbbkl = Ebıb Ebkl + Ebıbk Ebl + Ebıl Ebbk , which holds when the joint distribution of x and x
Gaussian. Thus, Eq. 8 is approximated byp̄i (1 − p̄i )(1 − 2p̄i )∆Jii = 0,22and therefore ∆Jii = 0 for all i. Hence, Eq. 7 can be approximated by−(1 − δij )((1 − 2p̄i )(1 − 2p̄j )Eij + p̄i p̄j (1 − p̄i )(1 − p̄j ) − Eij2 )(∆Jij + ∆Jji )−δij p̄i (1 − p̄i )(1 − 2p̄i )∆hi
= −δijX
X(δik (1 − 2p̄i )(1 − 2p̄l )Ebıl + 2(1 − δik )Ebıbk Ebl )∆Wkl − (1 − δij )(Ebıbk Ebl + Ebıl Ebbk )∆Wkl ,klkl

where δij is the Kronecker delta. For i 6= j, we have−((1 − 2p̄i )(1 − 2p̄j )Eij + p̄i p̄j (1 − p̄i )(1 − p̄j ) − Eij2 )(∆Jij + ∆Jji )= −X(Ebıbk Ebl + Ebıl Ebbk )∆Wkl .
klThus, ∆Jij is given byPl + Ebıl Eb
1
kl (Ebıbk Eb
bk )∆Wkl.
∆Jij = ∆Jji = (1 − δij )2 (1 − 2p̄i )(1 − 2p̄j )Eij + p̄i p̄j (1 − p̄i )(1 − p̄j ) − Eij2
Hence, we have
X ∂EijZ ∂Jmn
∂Eij
=
∂Wkl
∂Jmn ∂Wkl
mn≈= (1 − δij )




∂Jij
∂Jji
+
∂Wkl ∂Wkl((1 − 2p̄i )(1 − 2p̄j )Eij + p̄i p̄j (1 − p̄i )(1 − p̄j ) − Eij2 )(Ebıbk Ebl + Ebıl Ebbk )((1 − 2p̄i )(1 − 2p̄j )Eij + p̄i p̄j (1 − p̄i )(1 − p̄j ) − Eij2 )(1 − 2p̄i )(1 − 2p̄j )Eij + p̄i p̄j (1 − p̄i )(1 − p̄j ) − Eij2

= (1 − δij )(Ebıbk Ebl + Ebıl Ebbk ).
Assuming that only Wij affects the correlation Ebıj and that Ebıj is independent of Jkl and hk , weobtain∂EbıZj
∂Ebıj
≈
∂Wkl
∂Wkl
= Ebıj bkl − Ebıj Ebkl
≈ δik δjl ((1 − 2p̄i )(1 − 2p̄j )Ebıj + p̄i p̄j (1 − p̄i )(1 − p̄j ) − Ebı2j ).

23

Therefore, we findX ∂Cij∂log |C| =(C −1 )ji
∂Wkl
∂W
kl
ij
X ∂Eij(C −1 )ji
=
∂W
kl
ij
X
≈
(1 − δij )(Ebıbk Ebl + Ebıl Ebbk )(C −1 )jiij= VklCand∂log |D| =∂Wkl∂Dij −1D∂Wkl ji1≤i≤2N 1≤j≤2N
X ∂EijX ∂Eij(D−1 )ji +(D−1 )j+N i+N
=
∂W
∂W
kl
kl
ij
ij
X ∂Eib
X ∂Ebıj
+(D−1 )j+N i +(D−1 )j i+N
∂W
∂W
kl
kl
ij
ij
X
X
≈
(1 − δij )(Ebıbk Ebl + Ebıl Ebbk )(D−1 )ji +(1 − δij )(Ebıbk Ebl + Ebıl Ebbk )(D−1 )j+N i+N
X

ij

ij

+

X

+

X

δjk δil ((1 − 2p̄i )(1 − 2p̄j )Ebi + p̄i p̄j (1 − p̄i )(1 − p̄j ) − Eb2i )(D−1 )j+N iij

δik δjl ((1 − 2p̄i )(1 − 2p̄j )Ebıj + p̄i p̄j (1 − p̄i )(1 − p̄j ) − Ebı2j )(D−1 )j i+N

ij

=

X
(1 − δij )(Ebıbk Ebl + Ebıl Ebbk )((D−1 )ji + (D−1 )j+N i+N )ij

2)((D−1 )l k+N + (D−1 )k+N l )+((1 − 2p̄k )(1 − 2p̄l )Ebkl + p̄k p̄l (1 − p̄k )(1 − p̄l ) − Ebkl

= VklD .
Combining the above forms of∂∂log |C| andlog |D|, we find that the steepest gradient
∂Wkl
∂Wkl

24

Vkl of the approximate mutual information is given by∂b
I(X; X)
∂Wkl


1
∂
log |C| − log |D|
≈
∂Wkl
2
1
≈ VklC − VklD
2X
1=
(1 − δij )(Ebıbk Ebl + Ebıl Ebbk )(2(C −1 )ji − (D−1 )ji − (D−1 )j+N i+N )
2 ij

Vkl =

1
2
)((D−1 )l k+N + (D−1 )k+N l ).
− ((1 − 2p̄k )(1 − 2p̄l )Ebkl + p̄k p̄l (1 − p̄k )(1 − p̄l ) − Ebkl
2
To test the approximation of the above gradient, we compared the difference between ∆ log |C| =log |C 0 | − log |C| and0C
ij Vij ∆Wij , where C and C are the correlation matrices of x in the systems

Pwith the connection matrices Wij and Wij + ∆Wij , respectively. Fig. 11A shows that this approximation is quite good. The approximations ∆ log |D| = log |D0 | − log |D| ≈D
ij Vij ∆WijP

andb = I 0 (X; X)
b − I(X; X)b ≈ P Vij ∆Wij also give good estimates (Fig. 11B,C). Although∆I(X; X)ij
we set pmax = 1 for the system depicted in Fig. 11,P

ij Vij ∆Wij is a good index for the differenceb in the case pmax < 1.
of the mutual information ∆I(X; X),

B

Method for counting the repeated patterns and sequences

In Fig. 6, we present the number of repeated patterns and sequences before and after learning.
Defining repeated patterns as exact patterns that occur multiple times, we excluded incompletely
matched patterns from the definition of repeated patterns. This is because we wanted to simplify
the definition in order to make the result clear. Of course, we could define it such that patterns
with small differences can be regarded as a single repeated pattern. For example, if two patterns
with one mismatch or less, such as patterns a and b in Fig. 12A, are regarded as the same pattern,patterns b and c would also be regarded as the same pattern. Patterns a and c, however, cannot be25regarded as the same pattern, because the states of two of their neurons differ. Thus, in general,
even if some patterns, a and b, are considered to be the same pattern and some patterns, b and c,
are considered to be the same, patterns a and c may not be the same, according to this definition.
Thus, classifying two slightly different patterns into one repeated pattern makes the definition of
the repeated patterns less meaningful.
We defined a repeated sequence as an exact series of patterns that occurs more than once in a
block. A repeated sequence is thus composed of repeated patterns. Moreover, a repeated sequenceis composed of shorter repeated sequences. For example, each repeated sequence of length 4
contains 3 repeated sequences of length 2 (Fig. 12B). In general, a repeated sequence of length l1
contains l1 − l2 + 1 repeated sequences of length l2 < l1 . At first glance, it might seem that thisway of counting repeated sequences overestimates the number of occurrences of repeated sequences
and should be replaced by some more sophisticated method, such as a definition that does notcount the short sequences contained in a longer repeated sequence as a repeated sequence. Sucha method of counting, however, underestimates the number of repeated sequences. If a sequence
B of length 2 occurs three times, twice in a repeated sequence D of length 4 (B2 in D1 and B3
in D2 of Fig. 12C) and once outside of longer sequences (B1 in Fig. 12C), this modified way ofcounting fails to count the sequence B1 as an occurrence of the repeated sequence of length 2,even though this sequence is indeed repeated. To avoid this kind of failure, we counted sequences
as repeated sequences even when they were contained in longer repeated sequences. Thus, eachof the sequences A, B, C, and D occurs twice in Fig 12B, and the sequences B and D occur threetimes and twice, respectively, in Fig 12C.

Acknowledgments
This work was supported by Grants-in-Aid from the Ministry of Education, Science, Sports,

26

and Culture of Japan: Grant numbers 16200025, 17022020, 17650100, 18019019, 18047014, and18300079.

References
Abbott, L. and Rohrkemper, R. (2007). A simple growth model constructs critical avalanche
networks. Prog. Brain Res., 165:13–19.
Abeles, M. (1991). Corticonics. Cambridge. Univ. Press, Cambridge.
Atick, J. J. (1992). Could information theory provide an ecological theory of sensory processing?
Network, 3:213–251.
Ay, N. (2002). Locality of global stochastic interaction in directed acyclic networks. NeuralComput., 14(12):2959–2980.
Bak, P., Tang, C., and Wiesenfeld, K. (1987). Self-organized criticality: An explanation of the 1/f
noise. Phys. Rev. Lett., 59(4):381–384.
Beggs, J. M. and Plenz, D. (2003). Neuronal avalanches in neocortical circuits. J. Neurosci.,
23(35):11167–11177.
Bell, A. J. and Sejnowski, T. J. (1995). An information-maximization approach to blind separation
and blind deconvolution. Neural Comput., 7(6):1129–1159.
Bell, A. J. and Sejnowski, T. J. (1997). The ‘independent components’ of natural scenes are edge
filters. Vision Res., 37(23):3327–3338.
Buonomano, D. (2005). A learning rule for the emergence of stable dynamics and timing inrecurrent networks. J. Neurophysiol., 94:2275–2283.
27

Cossart, R., Aronov, D., and Yuste, R. (2003). Attractor dynamics of network UP states in the
neocortex. Nature, 423(6937):283–288.
Cover, T. M. and Thomas, J. A. (2006). Elements of Information Theory. John Wiley & Sons,
Inc., 2nd edition.
Diesmann, M., Gewaltig, M. O., and Aertsen, A. (1999). Stable propagation of synchronous spiking
in cortical neural networks. Nature, 402(6761):529–533.
Dragoi, G. and Buzsáki, G. (2006). Temporal encoding of place sequences by hippocampal cell
assemblies. Neuron, 50:145–157.
Fujisawa, S., Matsuki, N., and Ikegaya, Y. (2006). Single neurons can induce phase transitions ofcortical recurrent networks with multiple internal States. Cereb. Cortex, 16(5):639–654.
Harris, T. E. (1989). The theory of branching processes. Dover, New York.
Hebb, D. O. (1949). The Organization of Behavior; a Neuropsychological Theory. Wiley, NewYork.
Hinton, G. E. and Sejnowski, T. J. (1983). Optimal perceptual inference. Proceedings of the IEEE
conference on Computer Vision and Pattern Recognition, pages 448–453.
Hubel, D. H. and Wiesel, T. N. (1959). Receptive fields of single neurones in the cat’s striate
cortex. J. Physiol., 148:574–591.
Ikegaya, Y., Aaron, G., Cossart, R., Aronov, D., Lampl, I., Ferster, D., and Yuste, R. (2004).
Synfire chains and cortical songs: temporal modules of cortical activity. Science, 304(5670):559–
564.
Lewicki, M. S. (2002). Efficient coding of natural sounds. Nat. Neurosci., 5(4):356–363.
28

Linsker, R. (1988). Self-organization in a perceptual network. Computer, 21(3):105–117.
Louie, K. and Wilson, M. A. (2001). Temporally structured replay of awake hippocampal ensemble
activity during rapid eye movement sleep. Neuron, 29(1):145–156.
Maass, W., Natschläger, T., and Markram, H. (2002). Real-time computing without stable states:
a new framework for neural computation based on perturbations. Neural. Comput., 14:2531–
2560.
Nadasdy, Z., Hirase, H., Czurko, A., Csicsvari, J., and Buzsaki, G. (1999). Replay and timecompression of recurring spike sequences in the hippocampus. J. Neurosci., 19(21):9497–9507.
Olshausen, B. A. and Field, D. J. (1996). Emergence of simple-cell receptive field properties bylearning a sparse code for natural images. Nature, 381(6583):607–609.
Olshausen, B. A. and Field, D. J. (1997). Sparse coding with an overcomplete basis set: a strategy
employed by V1? Vision Res., 37(23):3311–3325.
Pfister, J., Toyoizumi, T., Barber, D., and Gerstner, W. (2006). Optimal spike-timing-dependent
plasticity for precise action potential firing in supervised learning. Neural Comput., 18:1318–
1348.
Sakurai, Y. and Takahashi, S. (2006). Dynamic synchrony of firing in the monkey prefrontal cortex
during working-memory tasks. J. Neurosci., 26(40):10141–10153.
Shannon, C. E. (1948). A mathematical theory of communication. Bell System Technical Journal,
27:379–423, 623–656.
Skaggs, W. E. and McNaughton, B. L. (1996). Replay of neuronal firing sequences in rat hippocampus during sleep following spatial experience. Science, 271(5257):1870–1873.
29

Teramae, J.-N. and Fukai, T. (2007). Local cortical circuit model inferred from power-law distributed neuronal avalanches. J. Comput. Neurosci., 22(3):301–312.
Toyoizumi, T., Pfister, J., Aihara, K., and Gerstner, W. (2005). Generalized Bienenstock-CooperMunro rule for spiking neurons that maximizes information transmission. Proc. Natl. Acad. Sci.
U.S.A., 102:5239–5244.
Tsukada, M., Ishii, N., and Sato, R. (1975). Temporal pattern discrimination of impulse sequencesin the computer-simulated nerve cells. Biol. Cybern., 17:19–28.
Vogels, T. and Abbott, L. (2005). Signal propagation and logic gating in networks of integrateand-fire neurons. J. Neurosci., 25:10786–10795.
Yao, H., Shi, L., Han, F., Gao, H., and Dan, Y. (2007). Rapid learning in cortical coding of visual
scenes. Nat. Neurosci., 10(6):772–778.

30

Figure 1: Feedforward network and recurrent network. (A) The feedforward network consideredby Linsker (1988) and Bell and Sejnowski (1995). An input signal at time step t is converted
to an output signal at time step t by the output neurons without memory. Infomax maximizes
the mutual information of the input and output. Feedback from output to input (dashed line)
is not considered in Linsker (1988) and Bell and Sejnowski (1995). (B,C) The recurrent network
investigated here. Recurrent connections convey the output of the last time step, and the outputof the neurons at time t is used as the input to the neurons at time t + 1. Their firing states
at time step t − 1 are synchronously updated to time step t. This network is equivalent to the
network depicted in A with feedback if the output at t in A is regarded as the state at t and usedas the input at t + 1. This network hence has memory of the previous state, and the informationis retained in the network. Although ‘Input’ and ‘Output’ are not included in B and C, replacing
the output of some neurons with input signal provides this network with an external input, andthe output of the neurons can be regarded as the output of the network.

31

Figure 2: Basic ideas of recurrent infomax. (A1,A2) A less-optimized and an optimized feedforward
network. Mutual information is calculated using the probability distribution of input and output.
Optimizing the connection weight maximizes the information transmitted from the input to the
output units. In other words, infomax minimizes the information loss. Bell and Sejnowski (1997)
reported that infomax in a feedforward network whose input consists of information taken froma photograph of a natural scene produces output units with simple cell-like selectivity. (B1,B2)
Recurrent infomax. Maximizing the mutual information of the two successive steps in the network
improves information retention and reduces the information loss over time.

32

Figure 3: Simulation process. (A) In each step, firing states of neurons are synchronously updated
using Eq. 1. A block consists of 20,000-100,000 steps. (B) Wij was updated using Eq. 3 at the endof every two blocks.

33

Figure 4: Mutual information of two successive states increases through learning process. In
the learning process consisting of 500 blocks, we observed increasing mutual information of two
successive states in the system whose behavior is shown in Fig. 5.

34

Figure 5: Formation of the feedforward structure through an algorithm based on RI in the model
network with external input. (A) The original photograph (1024×1024) of a pine tree was converted
to a gray-scaled, high-pass filtered image. Image patches (12 × 12) randomly selected from thehigh-pass filtered image were used as the external inputs to the network at each time step. (B1,C1)
Initially, 432 neurons were connected according to a random weight matrix. Of these neurons 144
were on-input, 144 were off-input, and 144 were output neurons. Each of the 144 pixels in an image
patch was linked to a pair of an on- and an off-input neuron in such a manner that the on-input
and off-input neurons were set to 1 (fire) only when the corresponding pixels had a positive andnegative sign, respectively. Output neurons fired spontaneously according to Eq. 1. The weightmatrix before learning is shown in C1. Initially, the connection weight Wij was a random matrix.
(B2,C2) After learning, feedforward structure from input to output neurons appeared in the model
network. (D1,2) Averaging the image patches that evoked firings of the output neurons revealed
that the output neurons, which did not exhibit clear selectivity before learning, responded to theGabor-like stimulus after learning. Here the following parameter values were used: N = 432,p̄input ≈ 0.15, p̄output = 0.002, pmax = 0.95,  = 0.01, η = 20, and wlimit = 100. Each of the 500
learning blocks consisted of 60,000 steps.

35

Figure 6:36

Repeated spatial patterns and spatiotemporal sequences occurred frequently in the networkwith pmax = 0.95 after learning. (A) Mutual information of two states interleaved with n − 1 steps,
I (n) , in 5th, 1200th, and 2500th blocks. I (n) is a monotonically decreasing function because twostates separated by n steps share less information than two states separated by n − 1 steps. I (n)takes a larger value in the optimized network after 2500 learning blocks than I (1) before learning.
In the 1200th block, I (n) has become larger than at the 5th block but have not been optimized yet.
Because I (1) is an approximation of mutual information of two successive states, I (1) takes quitea large value after learning. Although I (1) largely deviates from the mutual information, it is agood index of the information retained in a recurrent network (see Appendix A). (B) We define
a repeated pattern as a spatial firing pattern that is identically repeated at different time steps.
The size of a pattern is defined as the number of neurons firing in the pattern. A sequence that
contains a particular set of patterns appearing repeatedly in the same temporal order is called a“repeated sequence.” The size of a repeated sequence is defined as the sum of the sizes of thepatterns contained in it. (C) The numbers of occurrences of the patterns and sequences repeatedin the latter half of a test block (50,000 steps) were compared before and after learning. In thishistogram, only the sequences with sizes larger than 5l, where l is the length of the sequence,
were counted. (D1,2) When the repeated patterns in the latter 50,000 steps were colored, it was
found that no pattern occurred more than once in this short raster plot before learning (D1).
By contrast, several patterns appeared multiple times in the raster plot after learning (D2). In
addition, repeated sequences were found only in the raster plot after learning (red stars and blue
diamonds). (E1,2) The initial Wij with random weights (E1) evolved into a matrix with relatively
few strong weights (E2) after learning. Here the following parameter values were used: N = 50,p̄ = 0.05, pmax = 0.95,  = 0.01, η = 0.2, and wlimit = 100. Each of the 2500 learning blocks
consisted of 100,000 steps.

37Figure 7: The network in Fig. 6 (Optimized) and shuffled networks. (A) Approximate mutual
informations of two successive states of all 100 shuffled networks are smaller than that of theoriginal network after learning in Fig. 6 (I (1) = 39.7). (B) The numbers of occurrences of repeatedsequences with length 3 of all 100 shuffled networks are smaller than that of the original network
after learning in Fig. 6 (1606 occurrences).

38

Figure 8: Spontaneous activity of the recurrent network with pmax = 0.5 and pmax = 0.25. (A1,2)
Individual bursts in the spontaneous activity before (A1) and after learning (A2) for the networkwith pmax = 0.5 are indicated by different colors. The bursts before learning were short andfrequently interrupted by steps without firing, whereas the bursts after learning had much longer
durations. (B1,2) The initial Wij with random weights evolved into a matrix with relatively few
strong weights. Most rows and columns contained two strong excitatory connections (black dots);that is, most neurons had two strong inputs and two strong outputs. (C) Frequency distribution
P (s) of the burst size plotted as a function of the size, s. The black line corresponds to a slope of−1.5. Here the following parameter values were used: N = 50, p̄i = 0.01,  = 0.01, η = 0.2, andwlimit = 100. Each of the 15,000 learning blocks consisted of 20,000 steps.

39Figure 9: A feedforward structure was embedded in the model network by the temporallystructured stimulation. (A1,2) In the learning blocks, the state of neuron 1 was set to 1 (fire)at random intervals ranging from 50 to 99 steps. The first time step, t, is indicated by the arrowin A1. At t + 2, the state of neuron 3 was set to 1, and at t + 6, the state of neuron 2 was set to
1. In the test block after learning, only neuron 1 was set to 1 at random intervals ranging from 50
to 99 steps (A2). External stimulations are indicated by red circles. (B1, 2) The network activity
in an early learning block (B1) and the test block (B2). The steps at which neuron 1 was set to
1 are indicated by arrows, and externally evoked firings of neurons 1, 2, and 3 are indicated by
red circles. Although the states of neurons 2 and 3 were not set from the outside during the test
block, neurons 2 and 3 fired spontaneously six and two steps, respectively, after neuron 1 fired(as indicated by orange circles). The sequence of firings embedded by learning was replayed after
the spontaneous firing of neuron 1 (double arrows). (C) The weight matrix of the network after
learning (top) and its schematic representation (bottom) indicate a feedforward structure which
underlies the firing sequence starting from neuron 1 and containing neurons 3 and 2. Here thefollowing parameter values were used: N = 50, p̄i = 0.02, pmax = 0.98,  = 0.01, η = 1, andwlimit = 1000. Each of the 4250 learning blocks consisted of 60,000 steps, and we performed one
test block after 4250 learning blocks.

40

b and its approximation log |C| − log |D|/2. The mutual inforFigure 10: Comparison of I(X; X)
b is obtained through the direct measurement of the joint distribution P (x, xb) inmation, I(X; X),a block with 2 × 107 steps. Here the following parameter values were used: N = 4, pmax = 1,p̄i = 0.4, and  = 1.0 × 10−5 . The initial values of Wij were drawn from a uniform distribution ofb and its approximation to
the interval [−0.5, 0.5]. We converted the mutual information I(X; X)
bits by multiplying by the factor 1/ log 2.

41

b (C) and their approximaFigure 11: Comparison of ∆ log |C| (A), ∆ log |D| (B), and ∆I(X; X)
tions. The initial values of Wij were drawn from uniform distribution on [−0.5, 0.5], and thedifferences, ∆Wij , were drawn from a uniform distribution on [−0.005, 0.005]. Wij is changed tob are the differences
Wij + ∆Wij at the beginning of the block 2. ∆ log |C|, ∆ log |D|, and ∆I(X; X)
between the values measured in blocks 1 and 2. Each block consists of 2 × 107 steps. The samerandom number series were used in these blocks. Here the following parameter values were used:
N = 4, pmax = 1, p̄i = 0.4, and  = 1.0 × 10−5 .

42Figure 12: Counting the number of occurrences of the repeated sequences. (A) Only identical
repeats are counted as a repeat. Although the pattern b is very similar to pattern a and pattern c(one mismatch), they are not counted as repeated patterns. (B) Long repeated sequences containshorter repeated sequences. The repeated sequence D of length 4 contains the repeated sequences
A, B, and C of length 2. Two repeated sequences of length 3 are also contained in the repeatedsequence D. (C) If the sequences contained in a longer repeated sequence are not counted as
repeated sequences, the number of occurrences of the repeated sequences are underestimated.

43

A

B
Input

C

OutputInput at t

Output at t

= State at t–1= State at t

State

Feedback

Feedback

Figure 1: Feedforward network and recurrent network. (A) The feedforward network consideredby Linsker (1988) and Bell and Sejnowski (1995). An input signal at time step t is converted
to an output signal at time step t by the output neurons without memory. Infomax maximizes
the mutual information of the input and output. Feedback from output to input (dashed line)
is not considered in Linsker (1988) and Bell and Sejnowski (1995). (B,C) The recurrent network
investigated here. Recurrent connections convey the output of the last time step, and the outputof the neurons at time t is used as the input to the neurons at time t + 1. Their ﬁring states
at time step t − 1 are synchronously updated to time step t. This network is equivalent to the
network depicted in A with feedback if the output at t in A is regarded as the state at t and usedas the input at t + 1. This network hence has memory of the previous state, and the informationis retained in the network. Although ‘Input’ and ‘Output’ are not included in B and C, replacing
the output of some neurons with input signal provides this network with an external input, andthe output of the neurons can be regarded as the output of the network.

31

A1 Input X
Block 1
xi

..

xi+1

Output Y

B1

State X
Block 1..

^
Next state X

yi

x(t)

x(t+1)

yi+1

x(t+1)

x(t+2)..

..

Mutual information I(X;Y)^
Mutual information I(X;X)

Learning
between blocks

A2
Block b+1
xbT+i

..

xbT+i+1

B2
Block b+1

..

ybT+i

x(bT+t)x(bT+t+1)

ybT+i+1

x(bT+t+1)

x(bT+t+2)..

..

Mutual information I(X;Y)

Mutual information I(X;X)^

Figure 2: Basic ideas of recurrent infomax. (A1,A2) A less-optimized and an optimized feedforward
network. Mutual information is calculated using the probability distribution of input and output.
Optimizing the connection weight maximizes the information transmitted from the input to the
output units. In other words, infomax minimizes the information loss. Bell and Sejnowski (1997)
reported that infomax in a feedforward network whose input consists of information taken froma photograph of a natural scene produces output units with simple cell-like selectivity. (B1,B2)
Recurrent infomax. Maximizing the mutual information of the two successive steps in the network
improves information retention and reduces the information loss over time.
32

A
Block b−1

Block bBlock b+1

20,000∼100,000 steps = 1 block

{
...

...

Eq. 1

B

Eq. 1

Eq. 1

...

...

Eq. 1

Eq. 1

Eq. 1

Eq. 1

Eq. 1

Eq. 1

Eq. 1

Simulation = 500∼15,000 blocks{...

Update Wij (Eq. 3)

Block bUpdate Wij (Eq. 3)

Block b+1

Update Wij (Eq. 3)...

Figure 3: Simulation process. (A) In each step, ﬁring states of neurons are synchronously updated
using Eq. 1. A block consists of 20,000-100,000 steps. (B) Wij was updated using Eq. 3 at the endof every two blocks.

33

12

I(1) (bits)

10
8
6
4
2
0
0

100

200
300
Block

400

500

Figure 4: Mutual information of two successive states increases through learning process. In
the learning process consisting of 500 blocks, we observed increasing mutual information of two
successive states in the system whose behavior is shown in Fig. 5.

34

Original image

0.4

On-input neuron

Output neuron

On-input
Postsynaptic neuron

Off-input neuron

0.2
Off-input0—0.2

Output—0.4
On-input

Off-input

OutputPresynaptic neuron

Learning

Learning

Learning

Gray-scale & high-pass

80
70

Postsynaptic neuron

On-input60
50
40Off-input

30
20
10

Output0—10On-input Off-input
OutputPresynaptic neuron

Figure 5: Formation of the feedforward structure through an algorithm based on RI in the model
network with external input. (A) The original photograph (1024×1024) of a pine tree was converted
to a gray-scaled, high-pass ﬁltered image. Image patches (12 × 12) randomly selected from thehigh-pass ﬁltered image were used as the external inputs to the network at each time step. (B1,C1)
Initially, 432 neurons were connected according to a random weight matrix. Of these neurons 144
were on-input, 144 were oﬀ-input, and 144 were output neurons. Each of the 144 pixels in an image
patch was linked to a pair of an on- and an oﬀ-input neuron in such a manner that the on-input
and oﬀ-input neurons were set to 1 (ﬁre) only when the corresponding pixels had a positive andnegative sign, respectively. Output neurons ﬁred spontaneously according to Eq. 1. The weightmatrix before learning is shown in C1. Initially, the connection weight Wij was a random matrix.
(B2,C2) After learning, feedforward structure from input to output neurons appeared in the model
network. (D1,2) Averaging the image patches that evoked ﬁrings of the output neurons revealed
that the output neurons, which did not exhibit clear selectivity before learning, responded to theGabor-like stimulus after learning. Here the following parameter values were used: N = 432,p̄input ≈ 0.15, p̄output = 0.002, pmax = 0.95,  = 0.01, η = 20, and wlimit = 100. Each of the 500
learning blocks consisted of 60,000 steps.
35

B

∗

Repeated patten (size=4)(bits)

30

(n )

{

{

Repeated pattern (size=3)

5th block
1200th block
After learning (2500th block)

40

C

∗

Repeated sequence (length=2, size=7)50

step t t+1

t’

t’+1

I

20
10

...

0
2

4

6

8

...

...

10000
1000
100
10
1

5827
4887

10 12 14 16 18 20n

Before Learning
After Learning
5814

12

6
0

0

Pa
Pa
Se
Se
qu
tte
qu
tte
rn
en
rn
en
ce
ce
siz
siz
e>
e=
len
len
4
3o
gth
gth
r4
=1
=3
0

E1

D1
10
Neuron

1606

20
30
40

Postsynaptic neuron

A

0.4
10
0.2
20
0

30

–0.2

40

–0.4

50

50
0

50

100

150

200

Time stepD2

250

10 20 30 40 50
Presynaptic neuron★ Sequence A♦ Sequence B

★

♦ ★♦

★

Postsynaptic neuron

E2
Neuron

10
20
30
40

100
80
60
40
20
0
–20
–40
–60
–80
–100

10
20
30
40
50

50
0

50

100

150
Time step

Figure 6:

36

200

250

10 20 30 40 50
Presynaptic neuron

Repeated spatial patterns and spatiotemporal sequences occurred frequently in the networkwith pmax = 0.95 after learning. (A) Mutual information of two states interleaved with n − 1 steps,
I (n) , in 5th, 1200th, and 2500th blocks. I (n) is a monotonically decreasing function because twostates separated by n steps share less information than two states separated by n − 1 steps. I (n)takes a larger value in the optimized network after 2500 learning blocks than I (1) before learning.
In the 1200th block, I (n) has become larger than at the 5th block but have not been optimized yet.
Because I (1) is an approximation of mutual information of two successive states, I (1) takes quitea large value after learning. Although I (1) largely deviates from the mutual information, it is agood index of the information retained in a recurrent network (see Appendix A). (B) We deﬁnea repeated pattern as a spatial ﬁring pattern that is identically repeated at diﬀerent time steps.
The size of a pattern is deﬁned as the number of neurons ﬁring in the pattern. A sequence that
contains a particular set of patterns appearing repeatedly in the same temporal order is called a“repeated sequence.” The size of a repeated sequence is deﬁned as the sum of the sizes of thepatterns contained in it. (C) The numbers of occurrences of the patterns and sequences repeatedin the latter half of a test block (50,000 steps) were compared before and after learning. In thishistogram, only the sequences with sizes larger than 5l, where l is the length of the sequence,
were counted. (D1,2) When the repeated patterns in the latter 50,000 steps were colored, it was
found that no pattern occurred more than once in this short raster plot before learning (D1).
By contrast, several patterns appeared multiple times in the raster plot after learning (D2). In
addition, repeated sequences were found only in the raster plot after learning (red stars and blue
diamonds). (E1,2) The initial Wij with random weights (E1) evolved into a matrix with relatively
few strong weights (E2) after learning. Here the following parameter values were used: N = 50,p̄ = 0.05, pmax = 0.95,  = 0.01, η = 0.2, and wlimit = 100. Each of the 2500 learning blocks
consisted of 100,000 steps.

37

A 0.25

B 0.2
Shuffled

Frequency

Frequency

0.2
0.15
Optimized

0.1
0.05
0
15

20

25

30

35

40

45

Shuffled

0.18
0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0

Optimized

0

I(1)

500
1000
1500
Number of sequences

2000

Figure 7: The network in Fig. 6 (Optimized) and shuﬄed networks. (A) Approximate mutual
informations of two successive states of all 100 shuﬄed networks are smaller than that of theoriginal network after learning in Fig. 6 (I (1) = 39.7). (B) The numbers of occurrences of repeatedsequences with length 3 of all 100 shuﬄed networks are smaller than that of the original network
after learning in Fig. 6 (1606 occurrences).

38

Postsynaptic neuron

20
30
40

20
30
40
50

50
0

100

200
300
Time step400

B2
Postsynaptic neuron

10
20
30
40
50
0

100

200
300
Time step400-1

10 20 30 40 50
Presynaptic neuron

500

A2
Neuron

10

500

0

80
60

20

40

30

10
-2
10
-3
10
-4
10
-5
10
-6
10
-7
10
-8
10

1

10

100

1000

s

100
10

Before Learning (pmax=0.5)
After Learning (pmax=0.5)
After learning (pmax=0.25)

C 100

20

40

0

50

–20
10 20 30 40 50
Presynaptic neuron

P(s)Neuron

10

0.5
0.4
0.3
0.2
0.1
0
–0.1
–0.2
–0.3
–0.4
–0.5

P(s)

B1

A1

10
-1
10
-2
10
10-3
10-4
10-5
10-6
10-7
-8
10

Before Learning (pmax=0.5)
After Learning (pmax=0.5)
After learning (pmax=0.25)

10 20 30 40 50 60 70 80 90100
sFigure 8: Spontaneous activity of the recurrent network with pmax = 0.5 and pmax = 0.25. (A1,2)
Individual bursts in the spontaneous activity before (A1) and after learning (A2) for the networkwith pmax = 0.5 are indicated by diﬀerent colors. The bursts before learning were short andfrequently interrupted by steps without ﬁring, whereas the bursts after learning had much longer
durations. (B1,2) The initial Wij with random weights evolved into a matrix with relatively few
strong weights. Most rows and columns contained two strong excitatory connections (black dots);that is, most neurons had two strong inputs and two strong outputs. (C) Frequency distribution
P (s) of the burst size plotted as a function of the size, s. The black line corresponds to a slope of−1.5. Here the following parameter values were used: N = 50, p̄i = 0.01,  = 0.01, η = 0.2, andwlimit = 100. Each of the 15,000 learning blocks consisted of 20,000 steps.

39

A1

B1

Learning blocks
50-99 steps

C
1000

1

2
3

...
42-91 steps

Postsynaptic neuron

t+1 t+2 t+3 t+4 t+5 t+6

1

Neuron

t

7

10
20
30
40

500

34
3
5

0

49
18
11

–500

28

50
0

50

100

150

2

–1000
1

Time step

B2

Test block50-99 steps
t
1
2
3

28

10

t+1 t+2 t+3 t+4 t+5 t+6...
42-91 steps

NeuronA2

7 34 3 5 49 18 11 28 2
Presynaptic neuron

11

2

20

18

30

1

34

49

40

5
7

50
0

50

100

150

3

Time step

Figure 9: A feedforward structure was embedded in the model network by the temporallystructured stimulation. (A1,2) In the learning blocks, the state of neuron 1 was set to 1 (ﬁre)at random intervals ranging from 50 to 99 steps. The ﬁrst time step, t, is indicated by the arrowin A1. At t + 2, the state of neuron 3 was set to 1, and at t + 6, the state of neuron 2 was set to
1. In the test block after learning, only neuron 1 was set to 1 at random intervals ranging from 50
to 99 steps (A2). External stimulations are indicated by red circles. (B1, 2) The network activity
in an early learning block (B1) and the test block (B2). The steps at which neuron 1 was set to
1 are indicated by arrows, and externally evoked ﬁrings of neurons 1, 2, and 3 are indicated by
red circles. Although the states of neurons 2 and 3 were not set from the outside during the test
block, neurons 2 and 3 ﬁred spontaneously six and two steps, respectively, after neuron 1 ﬁred
(as indicated by orange circles). The sequence of ﬁrings embedded by learning was replayed after
the spontaneous ﬁring of neuron 1 (double arrows). (C) The weight matrix of the network after
learning (top) and its schematic representation (bottom) indicate a feedforward structure which
underlies the ﬁring sequence starting from neuron 1 and containing neurons 3 and 2. Here thefollowing parameter values were used: N = 50, p̄i = 0.02, pmax = 0.98,  = 0.01, η = 1, andwlimit = 1000. Each of the 4250 learning blocks consisted of 60,000 steps, and we performed one
test block after 4250 learning blocks.

40

Value function log|C|-log|D|/2

0.08
0.07
0.06
0.05
0.04
0.03
0.03

0.04

0.05

0.06

0.07

0.08

^

Mutual Information I(X ;X ) (bits)

 and its approximation log |C| − log |D|/2. The mutual inforFigure 10: Comparison of I(X; X)
 is obtained through the direct measurement of the joint distribution P (x, x) inmation, I(X; X),a block with 2 × 107 steps. Here the following parameter values were used: N = 4, pmax = 1,p̄i = 0.4, and  = 1.0 × 10−5 . The initial values of Wij were drawn from a uniform distribution of and its approximation to
the interval [−0.5, 0.5]. We converted the mutual information I(X; X)
bits by multiplying by the factor 1/ log 2.

41

Δ log|C|

0.0000500.004

C 0.001

Δ log|D|

0.002

Approximation

B

0.0001

Approximation

Approximation

A

0

^

Δ I(X;X) (bits)

0-0.002-0.00005-0.00005

00.00005

0.0001-0.002

Simulation

0
Simulation

0.002

0.004

0

0.001

Simulation

 (C) and their approximaFigure 11: Comparison of Δ log |C| (A), Δ log |D| (B), and ΔI(X; X)
tions. The initial values of Wij were drawn from uniform distribution on [−0.5, 0.5], and thediﬀerences, ΔWij , were drawn from a uniform distribution on [−0.005, 0.005]. Wij is changed to are the diﬀerences
Wij + ΔWij at the beginning of the block 2. Δ log |C|, Δ log |D|, and ΔI(X; X)
between the values measured in blocks 1 and 2. Each block consists of 2 × 107 steps. The samerandom number series were used in these blocks. Here the following parameter values were used:
N = 4, pmax = 1, p̄i = 0.4, and  = 1.0 × 10−5 .

42

A

B

D

{
{

D

...

{

{

B1

A

...

... ...

D1
B2

...

D2
B3

{{
{{

C

{

...

A

{

...

c

B

{

...

b

B

{

a

C

{

C

...

...

...

Figure 12: Counting the number of occurrences of the repeated sequences. (A) Only identical
repeats are counted as a repeat. Although the pattern b is very similar to pattern a and pattern c(one mismatch), they are not counted as repeated patterns. (B) Long repeated sequences containshorter repeated sequences. The repeated sequence D of length 4 contains the repeated sequences
A, B, and C of length 2. Two repeated sequences of length 3 are also contained in the repeatedsequence D. (C) If the sequences contained in a longer repeated sequence are not counted as
repeated sequences, the number of occurrences of the repeated sequences are underestimated.

43