# Neuronal Ensemble Decoding Using a Dynamical Maximum Entropy Model

Authors: Duho Sin; Jinsoo Kim; Jee Hyun Choi; Sung-Phil Kim
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

As advances in neurotechnology allow us to access the ensemble activity of multiple neurons simultaneously, many neurophysiologic studies have investigated how to decode neuronal ensemble activity. Neuronal ensemble activity from different brain regions exhibits a variety of characteristics, requiring substantially different decoding approaches. Among various models, a maximum entropy decoder is known to exploit not only individual firing activity but also interactions between neurons, extracting information more accurately for the cases with persistent neuronal activity and/or low-frequency firing activity. However, it does not consider temporal changes in neuronal states and therefore would be susceptible to poor performance for nonstationary neuronal information processing. To address this issue, we develop a novel decoder that extends a maximum entropy decoder to take time-varying neural information into account. This decoder blends a dynamical system model of neural networks into the maximum entropy model to better suit for nonstationary circumstances. From two simulation studies, we demonstrate that the proposed dynamic maximum entropy decoder could cope well with time-varying information, which the conventional maximum entropy decoder could not achieve. The results suggest that the proposed decoder may be able to infer neural information more effectively as it exploits dynamical properties of underlying neural networks.


---

## Body

## 1. Introduction

Ensemble data derived from neuronal population activities have been subject to numerous decoding attempts [1]. Traditional methods mainly focused on neural information of a single neuron averaged over multiple trials, consequently suffering from intertrial variations and neglecting ensemble information stemming from interactions between multiple neurons. Thus, decoding methods that directly tackle neuronal ensemble activity as a whole are more desirable in those terms [2].To this end, an ensemble decoding method called the population vector (PV) model was proposed where the ensemble state was represented as a weighted sum of the preferred directions of individual neurons’ firing rates [3]. Originally aimed for the analysis of primary motor cortical activity, the PV model was also utilized for the analysis of data from various regions such as primary visual cortex [1, 4]. The idea of the vector representation of neuronal ensemble activity was further extended in a form of the optimal linear estimator (OLE) [5].Both PV and OLE, categorized as structural analysis methods for ensemble data, can be applied when covariate events lie in a single space such as the direction of movement. But in many real cases, where such a condition is not met, a Bayesian decoding method may provide a better decoding solution. Rather than simply merging the ensemble data structurally, a Bayesian method chooses the event with the maximum a posteriori probability. This recursive Bayesian decoding method has been proposed for neuronal ensemble decoding in various forms such as the point process [6], the Kalman filter, or the Particle filter [7, 8]. By considering the fact that neuronal ensemble data contain vital information regarding correlations between multiple neurons [9], different Bayesian decoding approaches utilize such correlations directly for probability distribution estimation [10] via the statistical dynamics analysis [11–13], based on the maximum entropy principle. Being capable of extracting information from the data with a temporal resolution of millisecond [12], the maximum entropy decoding approach is known to be robust even to the neural data with a low firing rate [11, 14].The maximum entropy decoding method directly draws upon the stability of neuronal states under certain condition and is well known for its simplicity and consistent performance. It is especially known to well represent collective behavior of neural networks. However, it does not take temporal patterns of neuronal ensemble activity into consideration. Therefore, we propose a novel maximum entropy decoder that can embrace information regarding temporal dynamics of neuronal ensemble activity and thus enhance decoding performance. Especially, our new model is expected to play a vital role in the decoding of neuronal ensemble signals in prefrontal regions of the brain where persistence firing activity with a low firing rate is often observed.The paper is organized as follows. First we illustrate the basic concept of maximum entropy decoding for neuronal ensemble data and the computational models to implement this concept. Then, we describe an extended Ising decoder that has been recently proposed as an effective maximum entropy decoding method for various neuronal representations. Next, we propose a new decoder that incorporates temporal dynamics of neuronal ensemble into the maximum entropy decoding model. Simulation studies illustrating advantages of using the proposed model follow. Finally, we discuss advantages and limitations of the proposed model and possible future research directions.


## 2. Methods

### 2.1. Maximum Entropy Decoding

Let multivariate neuronal ensemble dataR be observed in response to stimuli S. Our objective is to determine which stimulus s
∈
S was applied by decoding a given neuronal data r
∈
R. One possible method would be to find s
∈
S that maximizes the conditional probability of s given r:
(1)
s
^
=
arg
max
⁡
s
∈
S
⁡
P
(
s
∣
r
)
.Due to the difficulty of completely describing such a conditional probability, the Bayes Rule is applied to represent the conditional probability in terms of likelihood and prior,(2)
s
^
=
arg
max
⁡
s
∈
S
⁡
P
(
r
∣
s
)
P
(
s
)
P
(
r
)
.In most cases, the likelihood ofP
(
r
∣
s
) describing a generative relationship between a specific stimulus and its consequential response can be determined. Also, since P
(
r
) is irrelevant to the inference of s
∈
S, the equation above can be reduced as follows:
(3)
s
^
=
arg
max
⁡
s
∈
S
⁡
P
(
r
∣
s
)
P
(
s
)
.
When estimating P
(
r
∣
s
), not only individual neuronal activity in response to a specific stimulus but also correlations among neurons should carry important information. By utilizing a measure that maximizes entropy, which can be applied with the information regarding neuronal correlations, we obtain a Bayesian model called the Ising model [14] and a method of decoding neuronal ensemble activity based on this model is called the Ising decoder [12]. The Ising decoder can be briefly described as follows.Let us assume that the firing rate of theith neuron at a given time instant is represented in a binary form: r
i
=
+
1 when fired and r
i
=
-
1 when silent. With a correlation between a specific stimulus s and the firing rate (Hz) of the ith neuron, denoted as {
h
i
}, together with {
J
i
j
} as the Lagrange multiplier and ∑
r
∈
R
P
(
r
∣
s
)
=
1 as the constraint, the probability that yields the maximum entropy is given by [12, 14]
(4)
P
(
r
∣
s
)
=
1
∑
r
∈
R
e
-
H
(
r
∣
s
)
e
-
H
(
r
∣
s
)
.
Here Hamiltonian H is defined by
(5)
H
(
r
∣
s
)
=
-
1
2
∑
i
,
j
J
i
j
(
s
)
r
i
r
j
-
∑
i
h
i
(
s
)
r
i
.
The denominator in (4) will be henceforth referred to as a partition function denoted as Z.In order to solve (4) with the Ising model, the partition function Z regarding s for given {
J
i
j
,
h
i
} must be calculated in advance. A direct approach, however, would require 2
N calculations for N neurons, rendering it intractable with a large number of neurons. Instead, one can use an approximation solution via the Markov chain Monte Carlo (MCMC) method [15, 16] or mean field approximation [17, 18]. In particular, mean field approximation of Z for given {
J
i
j
,
h
i
} can be expressed as [12]
(6)
log
⁡
⁡
Z
=
∑
i
log
⁡
⁡
(
2
cosh
⁡
⁡
(
h
i
+
W
i
)
)
-
∑
i
W
i
〈
r
i
〉
+
∑
i
<
j
J
i
j
〈
r
i
〉
〈
r
j
〉
.
Here, 〈
r
j
〉 indicates the mean firing rate of the jth neuron and W
i
=
∑
i
≠
j
J
i
j
〈
r
j
〉.The next problem would be to find{
J
i
j
,
h
i
} that minimizes the Hamiltonian H in order to obtain a maximum entropy distribution. The Thouless-Anderson-Palmer method can be utilized to approximate computation [11, 12, 18, 19]. Using this method we can find J
i
j that satisfies
(7)
(
C
-
1
)
i
j
=
-
J
i
j
-
2
J
i
j
2
〈
r
i
〉
〈
r
j
〉
,
where C
i
j
=
〈
r
i
r
j
〉
-
〈
r
i
〉
〈
r
j
〉. Then, h
i can be computed by
(8)
h
i
=
tanh
-
1
⁡
〈
r
i
〉
-
∑
i
J
i
j
〈
r
j
〉
+
〈
r
i
〉
∑
j
J
i
j
2
(
1
-
〈
r
j
〉
2
)
.


### 2.2. Extended Ising Decoder

One of the biggest advantages of using the Ising decoder is its temporal resolution of one-thousandth of a second [12]. But there may be occasions where a bigger bin size than 1 ms is required. In such cases, there would be multiple firings of a neuron within a single bin and thus resulting in information loss in the simple binary representation of firing activity (i.e., {
-
1
,
+
1
}). Hence, we should consider a system where a state of a neuron is represented in more diverse classes. Since the Ising model only takes binary states into consideration, a model with an extended flexibility is required; one of such models can be found in the Potts model.The Potts model is a generalized version of the Ising model that embraces a state space of∑
=
{
1,2
,
…
,
q
}
V, where an integer q
≥
2. Hamiltonian H for an arbitrary σ
∈
∑ is defined as
(9)
H
(
σ
)
=
-
J
∑
(
i
,
j
)
∈
E
δ
(
σ
i
,
σ
j
)
-
h
∑
i
∈
V
σ
i
.
Here, δ represents the kronecker delta function and σ
i is the firing state of a neuron i in the neuronal ensemble V. The subset of V, denoted as E, represents a set of neurons that are correlated with others. The Ising model is a special case of the Potts model with q
=
2, σ
i, σ
j
∈
{
-
1
,
+
1
}, and δ
(
σ
i
,
σ
j
)
=
(
1
/
2
)
(
1
+
σ
i
σ
j
) [20]. We can generalize the equation even further to reflect differing levels of interaction depending on the class of state such as
(10)
H
(
σ
)
=
-
1
2
∑
i
,
j
∈
V
∑
k
,
l
=
1
q
J
i
j
k
l
m
σ
i
,
k
m
σ
j
,
l
-
∑
i
∈
V
h
i
σ
i
,
where m
σ
i
,
r
=
q
δ
(
σ
i
r
)
-
1 [21].Representation of a neuronal state as a variable for the Potts model can be realized in multiple approaches: counting the total number of firings within a bin or measuring the ratio of a current firing rate with respect to a maximum firing rate for each neuron. The state space size would be equal to the maximum number of firings within a bin for the former, whereas the latter would yield an arbitrary size depending on a way of quantizing the firing rates.In the same line with the Ising model, the Potts model requires approximation of the partition function. Unfortunately, it is not easy to provide a closed-form solution for the partition function in the Potts model. Instead, we utilize the MCMC method as follows [15, 16, 22].Lets represent a Hamiltonian parameter pair (
J
,
h
) and let H
(
σ
;
s
) be H
(
σ
) conditioned on s; we would like to approximate Z
(
s
) via the MCMC method. Suppose that we know the information of Z
(
s
0
) for s
0
=
(
0,0
) and we want to estimate Z
(
s
1
) from it. If we define H
(
σ
;
s
)
=
H
(
σ
;
s
0
+
λ
(
s
1
-
s
0
)
), Z
(
s
λ
)
=
Z
(
s
0
+
λ
(
s
1
-
s
0
)
), then
(11)
log
⁡
⁡
Z
(
s
1
)
=
log
⁡
⁡
Z
(
s
0
)
+
∫
0
1
d
d
λ
log
⁡
⁡
Z
(
s
λ
)
d
λ
.
Taking a closer look inside the integral term:
(12)
d
d
λ
log
⁡
Z
(
s
λ
)
=
1
Z
(
s
λ
)
d
d
λ
Z
(
s
λ
)
=
1
Z
(
s
λ
)
∑
σ
∈
Σ
d
d
λ
e
-
H
(
σ
;
s
λ
)
=
1
Z
(
s
λ
)
∑
σ
∈
Σ
e
-
H
(
σ
;
s
λ
)
d
d
λ
(
-
H
(
σ
;
s
λ
)
)
=
E
(
d
d
λ
(
-
H
(
σ
;
s
λ
)
)
)
.
We can apply the mean of Markov chain generated by a Gibbs sampler to obtain the expected value in (12). For a sufficiently large K, the integral is approximated by the trapezoidal rule as
(13)
∫
0
1
d
d
λ
log
⁡
⁡
Z
(
s
λ
)
d
λ
≈
∑
n
=
0
K
n
K
+
1
E
(
d
d
λ
(
-
H
(
σ
;
s
n
/
K
)
)
)
.
Therefore,
(14)
Z
(
s
1
)
≈
Z
(
s
0
)
exp
⁡
⁡
(
∑
n
=
0
K
n
K
+
1
E
(
d
d
λ
(
-
H
(
σ
;
s
n
/
K
)
)
)
)
.The extended Ising decoder based on the Potts model will be henceforth referred to as the conventional maximum entropy decoder.


### 2.3. Dynamical Maximum Entropy Decoder

One of the neural network models exploiting statistical dynamics is the Hopfield network [23, 24]. It shows a distinct difference when compared to Perceptron [25, 26] models. Perceptron is aimed at physical systems and does not consider abstract population characteristics where the Hopfield model does. Moreover, while Perceptron is focused on synchronized systems, the Hopfield network can also be applied to more general dynamical systems.In the Hopfield network, a state of each neuron is represented asσ
∈
{
0,1
}, determined with a threshold θ and a synaptic connection strength J:
(15)
σ
i
=
{
1
,
if
∑
j
J
i
j
σ
j
+
h
i
≥
θ
i
,
0
,
if
∑
j
J
i
j
σ
j
+
h
i
<
θ
i
,
where i and j are index neurons and h represents an external input. The synaptic strength J is updated in a specific time range s as follows:
(16)
J
i
j
=
∑
s
(
2
σ
i
s
-
1
)
(
2
σ
j
s
-
1
)
.
In case of i
=
j, J
i
j
=
0.The most important feature is that if we assume symmetry such thatJ
i
j
=
J
j
i, the Lyapunov function exists for this neural network. The Lyapunov function, also referred to as an energy function, is given by
(17)
L
(
σ
)
=
-
1
2
∑
i
,
j
J
i
j
σ
i
σ
j
-
∑
i
h
i
σ
i
.The existence of the Lyapunov function is important because it analytically presents how the state of a system is likely to change. In such cases, the Lyapunov function is identical to Hamiltonian of the Ising model.The Hopfield network has been generalized based on the Potts model [21]. In this case, along with various degrees of interactions depending on diverse neural states, the Lyapunov function becomes
(18)
L
(
σ
)
=
-
1
2
∑
i
,
j
∈
V
∑
k
,
l
=
1
q
J
i
j
k
l
m
σ
i
,
k
m
σ
j
,
l
-
∑
i
∈
V
h
i
σ
i
,
where m
σ
i
,
r
=
q
δ
(
σ
i
r
)
-
1. Note that Hamiltonian in the Potts model can be obtained from this Lyapunov function.Now, we describe how to build a dynamical maximum entropy decoder. First, let us consider some basic aspects of a dynamical system. Let a dynamical system receiving an external inputs be described by
(19)
f
s
:
Σ
⟶
Σ
,
with a set of all possible neuronal states as Σ. We then make two basic assumptions on this dynamical system as follows.Assumption 1.
The dynamical system is assumed to approach a stable state over time with no input. Then, for everys
∈
S, σ
∈
Σ, and time t
∈
Z
+,
(20)
L
(
f
s
t
(
σ
)
)
≤
L
(
σ
)
.Assumption 2.
LetΩ
s
f
⊆
Σ be a pseudo-periodic set in response to an external input s. Then, we assume that
(21)
|
Ω
s
f
|
<
∞
,
f
s
n
(
σ
)
∈
Ω
s
f
.
Here n is an integer and can be sufficiently small. We also assume that an external stimulus s is static during its stimulation period.Next, we consider two cases when there is sufficient time for the system to be stabilized (Case 1) and when there is not (Case 2). For both cases, we need to estimateΩ
s
f with the assumptions above, which allow us to conjecture that any σ
1 and σ
2 are close to each other if σ
1
=
w
=
σ
2 for w
∈
Ω
s
f. Here w represents a pseudo-periodic discretized set belonging to Ω
s
f.


#### 2.3.1. Case 1

In this case, we first define a probability that a stateσ belongs to certain w:
(22)
P
(
σ
=
w
∈
Ω
s
f
)
.
Using Hamiltonian based on the Lyapunov function, we can compute a probability that σ belongs to w as
(23)
P
(
σ
=
w
)
=
P
w
(
σ
)
∑
w
′
P
w
′
(
σ
)
.The next step is to define a statistical state transition function. We represent every possible pair of states given a change of input froms to s
′ as Σ
~; that is, 〈
σ
t
,
σ
t
+
1
〉
=
σ
~
=
Σ
~. Then, the transition function is defined as
(24)
Φ
w
,
w
′
(
w
,
w
′
)
=
1
|
Σ
~
|
∑
σ
~
∈
Σ
~
P
(
σ
t
=
w
)
P
(
σ
t
+
1
=
w
′
)
.
Decoding an input S
t
+
1 with the knowledge of S
t from neuronal ensemble data is stated as the following problem:
(25)
S
t
+
1
^
=
arg
max
⁡
s
∈
S
⁡
P
(
S
t
+
1
=
s
∣
σ
t
,
σ
t
+
1
,
S
t
)
.
This problem can be solved using the transition function and probabilities described above:
(26)
P
(
S
t
+
1
=
s
∣
σ
t
,
σ
t
+
1
,
S
t
)
=
∑
w
t
+
1
∈
Ω
s
f
P
(
σ
t
+
1
=
w
t
+
1
)
×
(
∑
w
t
∈
Ω
s
f
P
(
σ
t
=
w
t
)
Φ
S
t
+
1
,
s
(
w
t
,
w
t
+
1
)
)
.Note that we can also use a simplified version of this formula given by(27)
P
(
S
t
+
1
=
s
∣
σ
t
,
σ
t
+
1
)
=
∑
w
t
+
1
,
w
t
∈
Ω
s
f
P
(
σ
t
=
w
t
)
P
(
σ
t
+
1
=
w
t
+
1
)
Φ
s
(
w
t
,
w
t
+
1
)
.


#### 2.3.2. Case 2

Let∑ be a set of all possible neuronal ensemble states. Suppose there is a function on ∑ such as
(28)
∅
:
∑
×
S
⟶
∑
.
∅ can be seen as the restrictive function f
|
∑
 of a dynamical system f and also as a reconstruction of f with respect to ∑. What this entails is the fact that f
(
σ
)
¯
=
〈
f
1
(
σ
)
,
f
2
(
σ
)
,
…
,
f
n
+
1
(
σ
)
〉 can be found if s
¯
∈
∏
i
=
0
n
S
i is given for all σ
∈
∑. Since solving the problem with every possible σ
∈
∑ and s
∈
S is difficult, we use a technique that limits the length of the sequence to N.Rather than estimatingΩ
s
f as defined, we estimate the neural trajectory of a specific σ. For a small enough N, we substitute ω for s
¯
∈
∏
i
=
0
n
S
i, therefore acquiring the probability P
(
σ
=
ω
) from (23). Now let Ω
=
{
ω
}, a temporary set that retains its elements but not further down the line. With this set, we redefine the function ∅ as
(29)
∅
:
Ω
×
S
⟶
Ω
and also in terms of Φ as
(30)
∅
(
Ω
,
s
∈
S
)
=
Φ
s
∈
S
(
Ω
,
Ω
)
∈
[
0,1
]
,
such that the condition of ∑
ω
,
ω
′
∈
Ω
Φ
s
(
ω
,
ω
′
)
=
1 is satisfied. Hence, now our goal is to find a dynamical function Φ over time. The initial value of Φ can be assumed to follow the uniform distribution.As time elapses,σ
¯
=
〈
σ
1
,
σ
2
,
…
,
σ
N
〉 and s
¯
=
ω turns into σ
′
¯
=
〈
σ
2
,
…
,
σ
N
,
σ
N
+
1
〉 and s
′
¯
=
ω
′. Then, we update Φ with the result of (26) as follows:
(31)
d
d
t
Φ
s
(
ω
,
ω
′
)
=
η
1
(
val
-
Φ
s
(
ω
,
ω
′
)
)
.
The “val” in above equation can be regarded as either 1 or 0. Once Φ
s is properly normalized, we accept and normalize the transformation with a probability η
2 if the updated P
(
σ
′
=
ω
′
) enhances the performance of (26). η
1 and η
2 introduced here are time-varying and converge to 0 when t
→
∞.This method is partially in line with the simulated annealing technique. Hence, it is expected to approach to a right solution but, on the other hand, can consume longer computational time. The main characteristic of this method is thatω is not limited only to stable states. If not enough time is given to reach a stable state after each external input, it might be better to consider the temporary states generated during the process. The proposed method realizes such temporary states. Even the trajectory of state transition as a result of change in external input can be represented in ω. Although each ω does not necessarily represent the stable state of a neural network in such cases, it can be advantageous for the decoding system that involves temporal dynamics since it can represent the most likely state in response to the change.


### 2.4. Simulation Procedure

A synthetic neuronal ensemble data set was generated in order to test the new decoding method for both cases, “Case 1” and “Case 2,” discussed above.


#### 2.4.1. Case 1

A model capable of representing persistence activity of neurons is generally complex [27]; thus, for the sake of the simplicity of experiment, we generated neuronal ensemble data with fixed firing rates of sixteen neurons for each external stimulus, based on Poisson process. Stimuli were limited to two classes, S
1 and S
2, where two stable states existed for each. In short, for each stimulus, S
=
{
S
1
,
S
2
}, stable states Ω
S
1
f
=
{
w
1,1
,
w
1,2
} and Ω
S
2
f
=
{
w
2,1
,
w
2,2
} exist (see Figure 2(a)).While four out of sixteen neurons represented the actual state, the rest were set to have random mean firing rates between 0 Hz to 10 Hz to represent irrelevant activity in the ensemble data. Those four state-representing neurons were divided into two pairs with neurons that were assigned the same mean firing rate. Thus, each state could be expressed as the following pairs:w
1,1
∶
=
(
10,10,0
,
0
) Hz, w
1,2
∶
=
(
0,0
,
10,10
) Hz, w
2,1
∶
=
(
0,0
,
0,0
) Hz, and w
2,2
∶
=
(
10,10,10,10
) Hz.


#### 2.4.2. Case 2

In order to generate neuronal ensemble data for Case 2, we utilized the orientation tuning model [28] that represents a biological neural network model for directional coding [29, 30]. Similar to an excitatory-inhibitory neural network model, the model contains an ensemble of neurons that respond to the angular direction θ and delivers excitatory or inhibitory signals among neurons depending on the directional similarity. In our simulation, we defined a signal passed on from one population of direction θ to another population of θ
′ as
(32)
J
θ
′
θ
=
-
J
0
+
J
2
cos
⁡
⁡
(
θ
-
θ
′
)
.
Here, J
0 and J
2 as a balancing parameter were set as J
0
=
0.14 and J
2
=
0.5 for our simulation. The actual input that the neurons tuned to θ
′ received from an external input with an angle θ was given by
(33)
h
θ
′
θ
ext
=
c
(
1
-
ϵ
+
ϵ
cos
⁡
⁡
(
θ
-
θ
′
)
)
,
where ϵ and c were parameters representing the degree of input focus and overall intensity of the input, respectively. The range of ϵ was given as [
0
,
1
/
2
]; ϵ
=
0 indicated that every neuron received uniform input regardless of their inherent direction, and ϵ
=
1
/
2 indicated that less input was given if θ
′ was more distant from θ. Hence, a condition |
θ
′
-
θ
|
=
π left neurons of θ
′ unaffected. For this simulation, ϵ
=
0.5 and c
=
0.4 were used (see Figure 1).Simulated data generated from the orientation tuning model. For each of 8 directional stimuli provided to each neuron (a), simulated neuronal responses are generated (b). Horizontal axis for each box represents time-flow (for 200 ms, from left to right), and vertical axis represents single neurons. Magnitude of stimulation and neuronal response is represented by a greyscale colormap.
(a)
(b)Ensemble state changes and its estimation. (a) Grey arrows represent change of statesw
i
,
j
∈
Ω
S
i
f for a specific stimulus S
k. (b) Estimated ensemble state transition of Φ
s
,
s
′
(
w
,
w
′
) for different external inputs. For each transition from S
i to S
j, the probability of state change from w
i
,
k to w
j
,
l is expressed in a greyscale mapping.
(a)
State changes resulting from external stimuli
(b)
An example of estimatedΦ
s
,
s
′
(
w
,
w
′
)A temporal change of a neural network was expressed as a firing rate based on the neural network model. The firing rateU
i of the ith neuron was
(34)
τ
i
d
U
i
d
t
=
-
U
i
+
S
i
(
∑
i
J
j
i
U
j
)
.
Here, τ
i was a time constant representing a synaptic delay between neurons, J
j
i was the synaptic strength between the ith and jth neurons, and S
i was a threshold function [31]. This part of the simulation utilized a total of 2
9
=
512 neurons as well as a synaptic delay constant, τ
=
10 ms. A threshold function for the firing rate was defined with the scale of 10
2 Hz:
(35)
S
(
U
)
=
{
0
,
if
U
<
0
,
U
,
if
0
≤
U
≤
4
,
4
,
if
U
>
4
.
The number of neurons for each of eight directions was equal (
2
9
-
3
=
64
) and the synaptic strength J and the external input h were identical to the setting of orientation tuning model of Case 1. The probability that neuronal connections between neurons exist was set below 0.5.


## 2.1. Maximum Entropy Decoding

Let multivariate neuronal ensemble dataR be observed in response to stimuli S. Our objective is to determine which stimulus s
∈
S was applied by decoding a given neuronal data r
∈
R. One possible method would be to find s
∈
S that maximizes the conditional probability of s given r:
(1)
s
^
=
arg
max
⁡
s
∈
S
⁡
P
(
s
∣
r
)
.Due to the difficulty of completely describing such a conditional probability, the Bayes Rule is applied to represent the conditional probability in terms of likelihood and prior,(2)
s
^
=
arg
max
⁡
s
∈
S
⁡
P
(
r
∣
s
)
P
(
s
)
P
(
r
)
.In most cases, the likelihood ofP
(
r
∣
s
) describing a generative relationship between a specific stimulus and its consequential response can be determined. Also, since P
(
r
) is irrelevant to the inference of s
∈
S, the equation above can be reduced as follows:
(3)
s
^
=
arg
max
⁡
s
∈
S
⁡
P
(
r
∣
s
)
P
(
s
)
.
When estimating P
(
r
∣
s
), not only individual neuronal activity in response to a specific stimulus but also correlations among neurons should carry important information. By utilizing a measure that maximizes entropy, which can be applied with the information regarding neuronal correlations, we obtain a Bayesian model called the Ising model [14] and a method of decoding neuronal ensemble activity based on this model is called the Ising decoder [12]. The Ising decoder can be briefly described as follows.Let us assume that the firing rate of theith neuron at a given time instant is represented in a binary form: r
i
=
+
1 when fired and r
i
=
-
1 when silent. With a correlation between a specific stimulus s and the firing rate (Hz) of the ith neuron, denoted as {
h
i
}, together with {
J
i
j
} as the Lagrange multiplier and ∑
r
∈
R
P
(
r
∣
s
)
=
1 as the constraint, the probability that yields the maximum entropy is given by [12, 14]
(4)
P
(
r
∣
s
)
=
1
∑
r
∈
R
e
-
H
(
r
∣
s
)
e
-
H
(
r
∣
s
)
.
Here Hamiltonian H is defined by
(5)
H
(
r
∣
s
)
=
-
1
2
∑
i
,
j
J
i
j
(
s
)
r
i
r
j
-
∑
i
h
i
(
s
)
r
i
.
The denominator in (4) will be henceforth referred to as a partition function denoted as Z.In order to solve (4) with the Ising model, the partition function Z regarding s for given {
J
i
j
,
h
i
} must be calculated in advance. A direct approach, however, would require 2
N calculations for N neurons, rendering it intractable with a large number of neurons. Instead, one can use an approximation solution via the Markov chain Monte Carlo (MCMC) method [15, 16] or mean field approximation [17, 18]. In particular, mean field approximation of Z for given {
J
i
j
,
h
i
} can be expressed as [12]
(6)
log
⁡
⁡
Z
=
∑
i
log
⁡
⁡
(
2
cosh
⁡
⁡
(
h
i
+
W
i
)
)
-
∑
i
W
i
〈
r
i
〉
+
∑
i
<
j
J
i
j
〈
r
i
〉
〈
r
j
〉
.
Here, 〈
r
j
〉 indicates the mean firing rate of the jth neuron and W
i
=
∑
i
≠
j
J
i
j
〈
r
j
〉.The next problem would be to find{
J
i
j
,
h
i
} that minimizes the Hamiltonian H in order to obtain a maximum entropy distribution. The Thouless-Anderson-Palmer method can be utilized to approximate computation [11, 12, 18, 19]. Using this method we can find J
i
j that satisfies
(7)
(
C
-
1
)
i
j
=
-
J
i
j
-
2
J
i
j
2
〈
r
i
〉
〈
r
j
〉
,
where C
i
j
=
〈
r
i
r
j
〉
-
〈
r
i
〉
〈
r
j
〉. Then, h
i can be computed by
(8)
h
i
=
tanh
-
1
⁡
〈
r
i
〉
-
∑
i
J
i
j
〈
r
j
〉
+
〈
r
i
〉
∑
j
J
i
j
2
(
1
-
〈
r
j
〉
2
)
.


## 2.2. Extended Ising Decoder

One of the biggest advantages of using the Ising decoder is its temporal resolution of one-thousandth of a second [12]. But there may be occasions where a bigger bin size than 1 ms is required. In such cases, there would be multiple firings of a neuron within a single bin and thus resulting in information loss in the simple binary representation of firing activity (i.e., {
-
1
,
+
1
}). Hence, we should consider a system where a state of a neuron is represented in more diverse classes. Since the Ising model only takes binary states into consideration, a model with an extended flexibility is required; one of such models can be found in the Potts model.The Potts model is a generalized version of the Ising model that embraces a state space of∑
=
{
1,2
,
…
,
q
}
V, where an integer q
≥
2. Hamiltonian H for an arbitrary σ
∈
∑ is defined as
(9)
H
(
σ
)
=
-
J
∑
(
i
,
j
)
∈
E
δ
(
σ
i
,
σ
j
)
-
h
∑
i
∈
V
σ
i
.
Here, δ represents the kronecker delta function and σ
i is the firing state of a neuron i in the neuronal ensemble V. The subset of V, denoted as E, represents a set of neurons that are correlated with others. The Ising model is a special case of the Potts model with q
=
2, σ
i, σ
j
∈
{
-
1
,
+
1
}, and δ
(
σ
i
,
σ
j
)
=
(
1
/
2
)
(
1
+
σ
i
σ
j
) [20]. We can generalize the equation even further to reflect differing levels of interaction depending on the class of state such as
(10)
H
(
σ
)
=
-
1
2
∑
i
,
j
∈
V
∑
k
,
l
=
1
q
J
i
j
k
l
m
σ
i
,
k
m
σ
j
,
l
-
∑
i
∈
V
h
i
σ
i
,
where m
σ
i
,
r
=
q
δ
(
σ
i
r
)
-
1 [21].Representation of a neuronal state as a variable for the Potts model can be realized in multiple approaches: counting the total number of firings within a bin or measuring the ratio of a current firing rate with respect to a maximum firing rate for each neuron. The state space size would be equal to the maximum number of firings within a bin for the former, whereas the latter would yield an arbitrary size depending on a way of quantizing the firing rates.In the same line with the Ising model, the Potts model requires approximation of the partition function. Unfortunately, it is not easy to provide a closed-form solution for the partition function in the Potts model. Instead, we utilize the MCMC method as follows [15, 16, 22].Lets represent a Hamiltonian parameter pair (
J
,
h
) and let H
(
σ
;
s
) be H
(
σ
) conditioned on s; we would like to approximate Z
(
s
) via the MCMC method. Suppose that we know the information of Z
(
s
0
) for s
0
=
(
0,0
) and we want to estimate Z
(
s
1
) from it. If we define H
(
σ
;
s
)
=
H
(
σ
;
s
0
+
λ
(
s
1
-
s
0
)
), Z
(
s
λ
)
=
Z
(
s
0
+
λ
(
s
1
-
s
0
)
), then
(11)
log
⁡
⁡
Z
(
s
1
)
=
log
⁡
⁡
Z
(
s
0
)
+
∫
0
1
d
d
λ
log
⁡
⁡
Z
(
s
λ
)
d
λ
.
Taking a closer look inside the integral term:
(12)
d
d
λ
log
⁡
Z
(
s
λ
)
=
1
Z
(
s
λ
)
d
d
λ
Z
(
s
λ
)
=
1
Z
(
s
λ
)
∑
σ
∈
Σ
d
d
λ
e
-
H
(
σ
;
s
λ
)
=
1
Z
(
s
λ
)
∑
σ
∈
Σ
e
-
H
(
σ
;
s
λ
)
d
d
λ
(
-
H
(
σ
;
s
λ
)
)
=
E
(
d
d
λ
(
-
H
(
σ
;
s
λ
)
)
)
.
We can apply the mean of Markov chain generated by a Gibbs sampler to obtain the expected value in (12). For a sufficiently large K, the integral is approximated by the trapezoidal rule as
(13)
∫
0
1
d
d
λ
log
⁡
⁡
Z
(
s
λ
)
d
λ
≈
∑
n
=
0
K
n
K
+
1
E
(
d
d
λ
(
-
H
(
σ
;
s
n
/
K
)
)
)
.
Therefore,
(14)
Z
(
s
1
)
≈
Z
(
s
0
)
exp
⁡
⁡
(
∑
n
=
0
K
n
K
+
1
E
(
d
d
λ
(
-
H
(
σ
;
s
n
/
K
)
)
)
)
.The extended Ising decoder based on the Potts model will be henceforth referred to as the conventional maximum entropy decoder.


## 2.3. Dynamical Maximum Entropy Decoder

One of the neural network models exploiting statistical dynamics is the Hopfield network [23, 24]. It shows a distinct difference when compared to Perceptron [25, 26] models. Perceptron is aimed at physical systems and does not consider abstract population characteristics where the Hopfield model does. Moreover, while Perceptron is focused on synchronized systems, the Hopfield network can also be applied to more general dynamical systems.In the Hopfield network, a state of each neuron is represented asσ
∈
{
0,1
}, determined with a threshold θ and a synaptic connection strength J:
(15)
σ
i
=
{
1
,
if
∑
j
J
i
j
σ
j
+
h
i
≥
θ
i
,
0
,
if
∑
j
J
i
j
σ
j
+
h
i
<
θ
i
,
where i and j are index neurons and h represents an external input. The synaptic strength J is updated in a specific time range s as follows:
(16)
J
i
j
=
∑
s
(
2
σ
i
s
-
1
)
(
2
σ
j
s
-
1
)
.
In case of i
=
j, J
i
j
=
0.The most important feature is that if we assume symmetry such thatJ
i
j
=
J
j
i, the Lyapunov function exists for this neural network. The Lyapunov function, also referred to as an energy function, is given by
(17)
L
(
σ
)
=
-
1
2
∑
i
,
j
J
i
j
σ
i
σ
j
-
∑
i
h
i
σ
i
.The existence of the Lyapunov function is important because it analytically presents how the state of a system is likely to change. In such cases, the Lyapunov function is identical to Hamiltonian of the Ising model.The Hopfield network has been generalized based on the Potts model [21]. In this case, along with various degrees of interactions depending on diverse neural states, the Lyapunov function becomes
(18)
L
(
σ
)
=
-
1
2
∑
i
,
j
∈
V
∑
k
,
l
=
1
q
J
i
j
k
l
m
σ
i
,
k
m
σ
j
,
l
-
∑
i
∈
V
h
i
σ
i
,
where m
σ
i
,
r
=
q
δ
(
σ
i
r
)
-
1. Note that Hamiltonian in the Potts model can be obtained from this Lyapunov function.Now, we describe how to build a dynamical maximum entropy decoder. First, let us consider some basic aspects of a dynamical system. Let a dynamical system receiving an external inputs be described by
(19)
f
s
:
Σ
⟶
Σ
,
with a set of all possible neuronal states as Σ. We then make two basic assumptions on this dynamical system as follows.Assumption 1.
The dynamical system is assumed to approach a stable state over time with no input. Then, for everys
∈
S, σ
∈
Σ, and time t
∈
Z
+,
(20)
L
(
f
s
t
(
σ
)
)
≤
L
(
σ
)
.Assumption 2.
LetΩ
s
f
⊆
Σ be a pseudo-periodic set in response to an external input s. Then, we assume that
(21)
|
Ω
s
f
|
<
∞
,
f
s
n
(
σ
)
∈
Ω
s
f
.
Here n is an integer and can be sufficiently small. We also assume that an external stimulus s is static during its stimulation period.Next, we consider two cases when there is sufficient time for the system to be stabilized (Case 1) and when there is not (Case 2). For both cases, we need to estimateΩ
s
f with the assumptions above, which allow us to conjecture that any σ
1 and σ
2 are close to each other if σ
1
=
w
=
σ
2 for w
∈
Ω
s
f. Here w represents a pseudo-periodic discretized set belonging to Ω
s
f.


### 2.3.1. Case 1

In this case, we first define a probability that a stateσ belongs to certain w:
(22)
P
(
σ
=
w
∈
Ω
s
f
)
.
Using Hamiltonian based on the Lyapunov function, we can compute a probability that σ belongs to w as
(23)
P
(
σ
=
w
)
=
P
w
(
σ
)
∑
w
′
P
w
′
(
σ
)
.The next step is to define a statistical state transition function. We represent every possible pair of states given a change of input froms to s
′ as Σ
~; that is, 〈
σ
t
,
σ
t
+
1
〉
=
σ
~
=
Σ
~. Then, the transition function is defined as
(24)
Φ
w
,
w
′
(
w
,
w
′
)
=
1
|
Σ
~
|
∑
σ
~
∈
Σ
~
P
(
σ
t
=
w
)
P
(
σ
t
+
1
=
w
′
)
.
Decoding an input S
t
+
1 with the knowledge of S
t from neuronal ensemble data is stated as the following problem:
(25)
S
t
+
1
^
=
arg
max
⁡
s
∈
S
⁡
P
(
S
t
+
1
=
s
∣
σ
t
,
σ
t
+
1
,
S
t
)
.
This problem can be solved using the transition function and probabilities described above:
(26)
P
(
S
t
+
1
=
s
∣
σ
t
,
σ
t
+
1
,
S
t
)
=
∑
w
t
+
1
∈
Ω
s
f
P
(
σ
t
+
1
=
w
t
+
1
)
×
(
∑
w
t
∈
Ω
s
f
P
(
σ
t
=
w
t
)
Φ
S
t
+
1
,
s
(
w
t
,
w
t
+
1
)
)
.Note that we can also use a simplified version of this formula given by(27)
P
(
S
t
+
1
=
s
∣
σ
t
,
σ
t
+
1
)
=
∑
w
t
+
1
,
w
t
∈
Ω
s
f
P
(
σ
t
=
w
t
)
P
(
σ
t
+
1
=
w
t
+
1
)
Φ
s
(
w
t
,
w
t
+
1
)
.


### 2.3.2. Case 2

Let∑ be a set of all possible neuronal ensemble states. Suppose there is a function on ∑ such as
(28)
∅
:
∑
×
S
⟶
∑
.
∅ can be seen as the restrictive function f
|
∑
 of a dynamical system f and also as a reconstruction of f with respect to ∑. What this entails is the fact that f
(
σ
)
¯
=
〈
f
1
(
σ
)
,
f
2
(
σ
)
,
…
,
f
n
+
1
(
σ
)
〉 can be found if s
¯
∈
∏
i
=
0
n
S
i is given for all σ
∈
∑. Since solving the problem with every possible σ
∈
∑ and s
∈
S is difficult, we use a technique that limits the length of the sequence to N.Rather than estimatingΩ
s
f as defined, we estimate the neural trajectory of a specific σ. For a small enough N, we substitute ω for s
¯
∈
∏
i
=
0
n
S
i, therefore acquiring the probability P
(
σ
=
ω
) from (23). Now let Ω
=
{
ω
}, a temporary set that retains its elements but not further down the line. With this set, we redefine the function ∅ as
(29)
∅
:
Ω
×
S
⟶
Ω
and also in terms of Φ as
(30)
∅
(
Ω
,
s
∈
S
)
=
Φ
s
∈
S
(
Ω
,
Ω
)
∈
[
0,1
]
,
such that the condition of ∑
ω
,
ω
′
∈
Ω
Φ
s
(
ω
,
ω
′
)
=
1 is satisfied. Hence, now our goal is to find a dynamical function Φ over time. The initial value of Φ can be assumed to follow the uniform distribution.As time elapses,σ
¯
=
〈
σ
1
,
σ
2
,
…
,
σ
N
〉 and s
¯
=
ω turns into σ
′
¯
=
〈
σ
2
,
…
,
σ
N
,
σ
N
+
1
〉 and s
′
¯
=
ω
′. Then, we update Φ with the result of (26) as follows:
(31)
d
d
t
Φ
s
(
ω
,
ω
′
)
=
η
1
(
val
-
Φ
s
(
ω
,
ω
′
)
)
.
The “val” in above equation can be regarded as either 1 or 0. Once Φ
s is properly normalized, we accept and normalize the transformation with a probability η
2 if the updated P
(
σ
′
=
ω
′
) enhances the performance of (26). η
1 and η
2 introduced here are time-varying and converge to 0 when t
→
∞.This method is partially in line with the simulated annealing technique. Hence, it is expected to approach to a right solution but, on the other hand, can consume longer computational time. The main characteristic of this method is thatω is not limited only to stable states. If not enough time is given to reach a stable state after each external input, it might be better to consider the temporary states generated during the process. The proposed method realizes such temporary states. Even the trajectory of state transition as a result of change in external input can be represented in ω. Although each ω does not necessarily represent the stable state of a neural network in such cases, it can be advantageous for the decoding system that involves temporal dynamics since it can represent the most likely state in response to the change.


## 2.3.1. Case 1

In this case, we first define a probability that a stateσ belongs to certain w:
(22)
P
(
σ
=
w
∈
Ω
s
f
)
.
Using Hamiltonian based on the Lyapunov function, we can compute a probability that σ belongs to w as
(23)
P
(
σ
=
w
)
=
P
w
(
σ
)
∑
w
′
P
w
′
(
σ
)
.The next step is to define a statistical state transition function. We represent every possible pair of states given a change of input froms to s
′ as Σ
~; that is, 〈
σ
t
,
σ
t
+
1
〉
=
σ
~
=
Σ
~. Then, the transition function is defined as
(24)
Φ
w
,
w
′
(
w
,
w
′
)
=
1
|
Σ
~
|
∑
σ
~
∈
Σ
~
P
(
σ
t
=
w
)
P
(
σ
t
+
1
=
w
′
)
.
Decoding an input S
t
+
1 with the knowledge of S
t from neuronal ensemble data is stated as the following problem:
(25)
S
t
+
1
^
=
arg
max
⁡
s
∈
S
⁡
P
(
S
t
+
1
=
s
∣
σ
t
,
σ
t
+
1
,
S
t
)
.
This problem can be solved using the transition function and probabilities described above:
(26)
P
(
S
t
+
1
=
s
∣
σ
t
,
σ
t
+
1
,
S
t
)
=
∑
w
t
+
1
∈
Ω
s
f
P
(
σ
t
+
1
=
w
t
+
1
)
×
(
∑
w
t
∈
Ω
s
f
P
(
σ
t
=
w
t
)
Φ
S
t
+
1
,
s
(
w
t
,
w
t
+
1
)
)
.Note that we can also use a simplified version of this formula given by(27)
P
(
S
t
+
1
=
s
∣
σ
t
,
σ
t
+
1
)
=
∑
w
t
+
1
,
w
t
∈
Ω
s
f
P
(
σ
t
=
w
t
)
P
(
σ
t
+
1
=
w
t
+
1
)
Φ
s
(
w
t
,
w
t
+
1
)
.


## 2.3.2. Case 2

Let∑ be a set of all possible neuronal ensemble states. Suppose there is a function on ∑ such as
(28)
∅
:
∑
×
S
⟶
∑
.
∅ can be seen as the restrictive function f
|
∑
 of a dynamical system f and also as a reconstruction of f with respect to ∑. What this entails is the fact that f
(
σ
)
¯
=
〈
f
1
(
σ
)
,
f
2
(
σ
)
,
…
,
f
n
+
1
(
σ
)
〉 can be found if s
¯
∈
∏
i
=
0
n
S
i is given for all σ
∈
∑. Since solving the problem with every possible σ
∈
∑ and s
∈
S is difficult, we use a technique that limits the length of the sequence to N.Rather than estimatingΩ
s
f as defined, we estimate the neural trajectory of a specific σ. For a small enough N, we substitute ω for s
¯
∈
∏
i
=
0
n
S
i, therefore acquiring the probability P
(
σ
=
ω
) from (23). Now let Ω
=
{
ω
}, a temporary set that retains its elements but not further down the line. With this set, we redefine the function ∅ as
(29)
∅
:
Ω
×
S
⟶
Ω
and also in terms of Φ as
(30)
∅
(
Ω
,
s
∈
S
)
=
Φ
s
∈
S
(
Ω
,
Ω
)
∈
[
0,1
]
,
such that the condition of ∑
ω
,
ω
′
∈
Ω
Φ
s
(
ω
,
ω
′
)
=
1 is satisfied. Hence, now our goal is to find a dynamical function Φ over time. The initial value of Φ can be assumed to follow the uniform distribution.As time elapses,σ
¯
=
〈
σ
1
,
σ
2
,
…
,
σ
N
〉 and s
¯
=
ω turns into σ
′
¯
=
〈
σ
2
,
…
,
σ
N
,
σ
N
+
1
〉 and s
′
¯
=
ω
′. Then, we update Φ with the result of (26) as follows:
(31)
d
d
t
Φ
s
(
ω
,
ω
′
)
=
η
1
(
val
-
Φ
s
(
ω
,
ω
′
)
)
.
The “val” in above equation can be regarded as either 1 or 0. Once Φ
s is properly normalized, we accept and normalize the transformation with a probability η
2 if the updated P
(
σ
′
=
ω
′
) enhances the performance of (26). η
1 and η
2 introduced here are time-varying and converge to 0 when t
→
∞.This method is partially in line with the simulated annealing technique. Hence, it is expected to approach to a right solution but, on the other hand, can consume longer computational time. The main characteristic of this method is thatω is not limited only to stable states. If not enough time is given to reach a stable state after each external input, it might be better to consider the temporary states generated during the process. The proposed method realizes such temporary states. Even the trajectory of state transition as a result of change in external input can be represented in ω. Although each ω does not necessarily represent the stable state of a neural network in such cases, it can be advantageous for the decoding system that involves temporal dynamics since it can represent the most likely state in response to the change.


## 2.4. Simulation Procedure

A synthetic neuronal ensemble data set was generated in order to test the new decoding method for both cases, “Case 1” and “Case 2,” discussed above.


### 2.4.1. Case 1

A model capable of representing persistence activity of neurons is generally complex [27]; thus, for the sake of the simplicity of experiment, we generated neuronal ensemble data with fixed firing rates of sixteen neurons for each external stimulus, based on Poisson process. Stimuli were limited to two classes, S
1 and S
2, where two stable states existed for each. In short, for each stimulus, S
=
{
S
1
,
S
2
}, stable states Ω
S
1
f
=
{
w
1,1
,
w
1,2
} and Ω
S
2
f
=
{
w
2,1
,
w
2,2
} exist (see Figure 2(a)).While four out of sixteen neurons represented the actual state, the rest were set to have random mean firing rates between 0 Hz to 10 Hz to represent irrelevant activity in the ensemble data. Those four state-representing neurons were divided into two pairs with neurons that were assigned the same mean firing rate. Thus, each state could be expressed as the following pairs:w
1,1
∶
=
(
10,10,0
,
0
) Hz, w
1,2
∶
=
(
0,0
,
10,10
) Hz, w
2,1
∶
=
(
0,0
,
0,0
) Hz, and w
2,2
∶
=
(
10,10,10,10
) Hz.


### 2.4.2. Case 2

In order to generate neuronal ensemble data for Case 2, we utilized the orientation tuning model [28] that represents a biological neural network model for directional coding [29, 30]. Similar to an excitatory-inhibitory neural network model, the model contains an ensemble of neurons that respond to the angular direction θ and delivers excitatory or inhibitory signals among neurons depending on the directional similarity. In our simulation, we defined a signal passed on from one population of direction θ to another population of θ
′ as
(32)
J
θ
′
θ
=
-
J
0
+
J
2
cos
⁡
⁡
(
θ
-
θ
′
)
.
Here, J
0 and J
2 as a balancing parameter were set as J
0
=
0.14 and J
2
=
0.5 for our simulation. The actual input that the neurons tuned to θ
′ received from an external input with an angle θ was given by
(33)
h
θ
′
θ
ext
=
c
(
1
-
ϵ
+
ϵ
cos
⁡
⁡
(
θ
-
θ
′
)
)
,
where ϵ and c were parameters representing the degree of input focus and overall intensity of the input, respectively. The range of ϵ was given as [
0
,
1
/
2
]; ϵ
=
0 indicated that every neuron received uniform input regardless of their inherent direction, and ϵ
=
1
/
2 indicated that less input was given if θ
′ was more distant from θ. Hence, a condition |
θ
′
-
θ
|
=
π left neurons of θ
′ unaffected. For this simulation, ϵ
=
0.5 and c
=
0.4 were used (see Figure 1).Simulated data generated from the orientation tuning model. For each of 8 directional stimuli provided to each neuron (a), simulated neuronal responses are generated (b). Horizontal axis for each box represents time-flow (for 200 ms, from left to right), and vertical axis represents single neurons. Magnitude of stimulation and neuronal response is represented by a greyscale colormap.
(a)
(b)Ensemble state changes and its estimation. (a) Grey arrows represent change of statesw
i
,
j
∈
Ω
S
i
f for a specific stimulus S
k. (b) Estimated ensemble state transition of Φ
s
,
s
′
(
w
,
w
′
) for different external inputs. For each transition from S
i to S
j, the probability of state change from w
i
,
k to w
j
,
l is expressed in a greyscale mapping.
(a)
State changes resulting from external stimuli
(b)
An example of estimatedΦ
s
,
s
′
(
w
,
w
′
)A temporal change of a neural network was expressed as a firing rate based on the neural network model. The firing rateU
i of the ith neuron was
(34)
τ
i
d
U
i
d
t
=
-
U
i
+
S
i
(
∑
i
J
j
i
U
j
)
.
Here, τ
i was a time constant representing a synaptic delay between neurons, J
j
i was the synaptic strength between the ith and jth neurons, and S
i was a threshold function [31]. This part of the simulation utilized a total of 2
9
=
512 neurons as well as a synaptic delay constant, τ
=
10 ms. A threshold function for the firing rate was defined with the scale of 10
2 Hz:
(35)
S
(
U
)
=
{
0
,
if
U
<
0
,
U
,
if
0
≤
U
≤
4
,
4
,
if
U
>
4
.
The number of neurons for each of eight directions was equal (
2
9
-
3
=
64
) and the synaptic strength J and the external input h were identical to the setting of orientation tuning model of Case 1. The probability that neuronal connections between neurons exist was set below 0.5.


## 2.4.1. Case 1

A model capable of representing persistence activity of neurons is generally complex [27]; thus, for the sake of the simplicity of experiment, we generated neuronal ensemble data with fixed firing rates of sixteen neurons for each external stimulus, based on Poisson process. Stimuli were limited to two classes, S
1 and S
2, where two stable states existed for each. In short, for each stimulus, S
=
{
S
1
,
S
2
}, stable states Ω
S
1
f
=
{
w
1,1
,
w
1,2
} and Ω
S
2
f
=
{
w
2,1
,
w
2,2
} exist (see Figure 2(a)).While four out of sixteen neurons represented the actual state, the rest were set to have random mean firing rates between 0 Hz to 10 Hz to represent irrelevant activity in the ensemble data. Those four state-representing neurons were divided into two pairs with neurons that were assigned the same mean firing rate. Thus, each state could be expressed as the following pairs:w
1,1
∶
=
(
10,10,0
,
0
) Hz, w
1,2
∶
=
(
0,0
,
10,10
) Hz, w
2,1
∶
=
(
0,0
,
0,0
) Hz, and w
2,2
∶
=
(
10,10,10,10
) Hz.


## 2.4.2. Case 2

In order to generate neuronal ensemble data for Case 2, we utilized the orientation tuning model [28] that represents a biological neural network model for directional coding [29, 30]. Similar to an excitatory-inhibitory neural network model, the model contains an ensemble of neurons that respond to the angular direction θ and delivers excitatory or inhibitory signals among neurons depending on the directional similarity. In our simulation, we defined a signal passed on from one population of direction θ to another population of θ
′ as
(32)
J
θ
′
θ
=
-
J
0
+
J
2
cos
⁡
⁡
(
θ
-
θ
′
)
.
Here, J
0 and J
2 as a balancing parameter were set as J
0
=
0.14 and J
2
=
0.5 for our simulation. The actual input that the neurons tuned to θ
′ received from an external input with an angle θ was given by
(33)
h
θ
′
θ
ext
=
c
(
1
-
ϵ
+
ϵ
cos
⁡
⁡
(
θ
-
θ
′
)
)
,
where ϵ and c were parameters representing the degree of input focus and overall intensity of the input, respectively. The range of ϵ was given as [
0
,
1
/
2
]; ϵ
=
0 indicated that every neuron received uniform input regardless of their inherent direction, and ϵ
=
1
/
2 indicated that less input was given if θ
′ was more distant from θ. Hence, a condition |
θ
′
-
θ
|
=
π left neurons of θ
′ unaffected. For this simulation, ϵ
=
0.5 and c
=
0.4 were used (see Figure 1).Simulated data generated from the orientation tuning model. For each of 8 directional stimuli provided to each neuron (a), simulated neuronal responses are generated (b). Horizontal axis for each box represents time-flow (for 200 ms, from left to right), and vertical axis represents single neurons. Magnitude of stimulation and neuronal response is represented by a greyscale colormap.
(a)
(b)Ensemble state changes and its estimation. (a) Grey arrows represent change of statesw
i
,
j
∈
Ω
S
i
f for a specific stimulus S
k. (b) Estimated ensemble state transition of Φ
s
,
s
′
(
w
,
w
′
) for different external inputs. For each transition from S
i to S
j, the probability of state change from w
i
,
k to w
j
,
l is expressed in a greyscale mapping.
(a)
State changes resulting from external stimuli
(b)
An example of estimatedΦ
s
,
s
′
(
w
,
w
′
)A temporal change of a neural network was expressed as a firing rate based on the neural network model. The firing rateU
i of the ith neuron was
(34)
τ
i
d
U
i
d
t
=
-
U
i
+
S
i
(
∑
i
J
j
i
U
j
)
.
Here, τ
i was a time constant representing a synaptic delay between neurons, J
j
i was the synaptic strength between the ith and jth neurons, and S
i was a threshold function [31]. This part of the simulation utilized a total of 2
9
=
512 neurons as well as a synaptic delay constant, τ
=
10 ms. A threshold function for the firing rate was defined with the scale of 10
2 Hz:
(35)
S
(
U
)
=
{
0
,
if
U
<
0
,
U
,
if
0
≤
U
≤
4
,
4
,
if
U
>
4
.
The number of neurons for each of eight directions was equal (
2
9
-
3
=
64
) and the synaptic strength J and the external input h were identical to the setting of orientation tuning model of Case 1. The probability that neuronal connections between neurons exist was set below 0.5.


## 3. Results and Discussion

Performance of the decoders was assessed by means of “decoding error rate” and “decoding uncertainty” for the given dataset. Error rate calculated as a ratio of incorrectly classified samples to the total number of samples tested served as a direct measure of decoder performance. Decoding uncertainty was represented with the entropy of the normalized distribution of Bayesian probabilities for each class generated by the decoder. Small amount of entropy indicates less uncertainty of decoding performance, analogous to smaller variance of estimates of a model. Using these performance measures, we compared our proposed decoder (dynamical maximum entropy decoder) with the conventional decoder (the extended Ising decoder).


### 3.1. Case 1

The dynamical maximum entropy decoder estimatedΦ
s
,
s
′
(
w
,
w
′
) as the first step of decoding. We examin ed whether this transition estimation was correctly represented by the decoder. Figure 2(b) demonstrates that the estimation results well represent state changes for a given external event (see Figure 2(a)). When there were more than two stable states for each event, the proposed method outperformed the conventional maximum entropy decoder in terms of both error rate and uncertainty (see Figure 3). The proposed decoder generated approximately 0.25 error rate and 0.12 bit of entropy which is a significant improvement over 0.45 error rate and 0.78 bit entropy by the conventional decoder. Superior performance of the dynamical maximum entropy decoder was based on its ability to alleviate signals from those neurons which were regarded as irrelevant to the event.(a) Decoding error rates and (b) uncertainty represented by entropy of decoding output using two decoding models from simulated neural data (see text). The maximum entropy is 1 bit in this case since there are two possible external events. In both categories, there is significant difference in error rates (t-test, P
<
0.00001) between the dynamical maximum entropy decoder (new) and the traditional maximum entropy decoder (traditional).
(a)
Decoding error
(b)
Decoding uncertainty


### 3.2. Case 2

The second study considered a case where the external stimulus changes before the system converges to a stable state. For each of 8 uniformly separated directions, the stimulus was changed after the first 200 ms and then lasted for another 200 ms period. The succeeding direction could either be different from or same as the former. A small number of neurons were randomly chosen out of 512 total neurons and were discretized within a bin size of 10 ms (see Figure4).Representation of sampled and discretized simulation data. 64 random neurons were sampled from virtual ensemble data with 512 neurons (a) and then discretized in a 10 ms resolution (100 Hz) (b). Stimulation direction was0.5
π for the first 120 ms and π for the next 240 ms.
(a)
Sampled simulation data
(b)
Discretized simulation dataFor these data with temporal changes of stimuli, the proposed decoder was evidently superior compared to the conventional decoder in terms of error rates with 0.07-, 0.07-, 0.01-, and 0.11-point advantage for respective scenarios (see Figure5(a)). However, decoding uncertainty was also significantly (difference range 0.3~0.8 bits; t-test, P
<
0.01) higher in the proposed decoder than the conventional one in all scenarios (Figure 5(b)). This issue will be further discussed below.A comparison of decoding error and uncertainty under various degrees of stimulus direction difference between the proposed (new) and the traditional maximum entropy decoders. (a) Both decoders demonstrate a similar trend of higher error rates for greater input difference. Except for the0.75
π condition, the proposed method yields significantly lower error rates (t-test, P
<
0.00001). (b) A sensible increase of uncertainty is shown for both methods as the input difference grows. The proposed method contains greater uncertainty (t-test, P
<
0.01) for all cases.
(a)
Decoding error
(b)
Decoding uncertainty


### 3.3. Discussion

Decoding neuronal ensemble activity recorded from brain signals provides insights on how neural networks process information. Decoding models based on computational models for biological neural networks can provide these insights to deepen our understanding of neural network mechanisms. However, approaches directly using neural network models are not easily realized because it is challenging to know every synaptic strength between neurons and to obtain full information of all the neurons in a network. Therefore, the present study proposes an alternative approach to leveraging dynamical and collective behavior of an ensemble of neurons in response to various stimuli. From simple assumptions on the model, we blend neural dynamics into a decoder with which we are able to inspect the functions and roles of a neural network.To the assessment of the performance of a decoder, the uncertainty of decoding outcome as well as the error rate is crucial. This uncertainty measures how robust the decision made by a decoder is regarding the information inferred from neuronal activity. Our proposed decoder exhibited a significantly low error rate and small uncertainty in the first case study (see Figure3). However, it generated significantly larger uncertainty than the conventional maximum entropy decoder in the second case study despite its lower error rate (Figure 5). We speculate that such larger uncertainty may be due to model complexity in our decoder. More complex decoding procedure in our model includes uniformly distributed transition probabilities that may in turn equalize prior probabilities of individual stimuli, thus increasing uncertainty of a decoded stimulus. However, it may not be solely the result of model structure because the proposed decoder could also reduce uncertainty in the first case. Consequently, we suspect that increased uncertainty in the second case may indicate a particular outcome resulting from specific data properties, which need further investigation in the future.Yet, we also recognize that there is plenty of room to improve in our model. In particular, more rigorous ways of obtaining information about neural dynamics may be necessary. For instance, our method to estimate stabilized states of a neuronal ensemble may be suboptimal to many different real data. Also, we need to apply our decoder to ensemble data from many different brain regions to generalize it further. Finally, continuous efforts to reduce computational loads in the proposed decoder will be required.


## 3.1. Case 1

The dynamical maximum entropy decoder estimatedΦ
s
,
s
′
(
w
,
w
′
) as the first step of decoding. We examin ed whether this transition estimation was correctly represented by the decoder. Figure 2(b) demonstrates that the estimation results well represent state changes for a given external event (see Figure 2(a)). When there were more than two stable states for each event, the proposed method outperformed the conventional maximum entropy decoder in terms of both error rate and uncertainty (see Figure 3). The proposed decoder generated approximately 0.25 error rate and 0.12 bit of entropy which is a significant improvement over 0.45 error rate and 0.78 bit entropy by the conventional decoder. Superior performance of the dynamical maximum entropy decoder was based on its ability to alleviate signals from those neurons which were regarded as irrelevant to the event.(a) Decoding error rates and (b) uncertainty represented by entropy of decoding output using two decoding models from simulated neural data (see text). The maximum entropy is 1 bit in this case since there are two possible external events. In both categories, there is significant difference in error rates (t-test, P
<
0.00001) between the dynamical maximum entropy decoder (new) and the traditional maximum entropy decoder (traditional).
(a)
Decoding error
(b)
Decoding uncertainty


## 3.2. Case 2

The second study considered a case where the external stimulus changes before the system converges to a stable state. For each of 8 uniformly separated directions, the stimulus was changed after the first 200 ms and then lasted for another 200 ms period. The succeeding direction could either be different from or same as the former. A small number of neurons were randomly chosen out of 512 total neurons and were discretized within a bin size of 10 ms (see Figure4).Representation of sampled and discretized simulation data. 64 random neurons were sampled from virtual ensemble data with 512 neurons (a) and then discretized in a 10 ms resolution (100 Hz) (b). Stimulation direction was0.5
π for the first 120 ms and π for the next 240 ms.
(a)
Sampled simulation data
(b)
Discretized simulation dataFor these data with temporal changes of stimuli, the proposed decoder was evidently superior compared to the conventional decoder in terms of error rates with 0.07-, 0.07-, 0.01-, and 0.11-point advantage for respective scenarios (see Figure5(a)). However, decoding uncertainty was also significantly (difference range 0.3~0.8 bits; t-test, P
<
0.01) higher in the proposed decoder than the conventional one in all scenarios (Figure 5(b)). This issue will be further discussed below.A comparison of decoding error and uncertainty under various degrees of stimulus direction difference between the proposed (new) and the traditional maximum entropy decoders. (a) Both decoders demonstrate a similar trend of higher error rates for greater input difference. Except for the0.75
π condition, the proposed method yields significantly lower error rates (t-test, P
<
0.00001). (b) A sensible increase of uncertainty is shown for both methods as the input difference grows. The proposed method contains greater uncertainty (t-test, P
<
0.01) for all cases.
(a)
Decoding error
(b)
Decoding uncertainty


## 3.3. Discussion

Decoding neuronal ensemble activity recorded from brain signals provides insights on how neural networks process information. Decoding models based on computational models for biological neural networks can provide these insights to deepen our understanding of neural network mechanisms. However, approaches directly using neural network models are not easily realized because it is challenging to know every synaptic strength between neurons and to obtain full information of all the neurons in a network. Therefore, the present study proposes an alternative approach to leveraging dynamical and collective behavior of an ensemble of neurons in response to various stimuli. From simple assumptions on the model, we blend neural dynamics into a decoder with which we are able to inspect the functions and roles of a neural network.To the assessment of the performance of a decoder, the uncertainty of decoding outcome as well as the error rate is crucial. This uncertainty measures how robust the decision made by a decoder is regarding the information inferred from neuronal activity. Our proposed decoder exhibited a significantly low error rate and small uncertainty in the first case study (see Figure3). However, it generated significantly larger uncertainty than the conventional maximum entropy decoder in the second case study despite its lower error rate (Figure 5). We speculate that such larger uncertainty may be due to model complexity in our decoder. More complex decoding procedure in our model includes uniformly distributed transition probabilities that may in turn equalize prior probabilities of individual stimuli, thus increasing uncertainty of a decoded stimulus. However, it may not be solely the result of model structure because the proposed decoder could also reduce uncertainty in the first case. Consequently, we suspect that increased uncertainty in the second case may indicate a particular outcome resulting from specific data properties, which need further investigation in the future.Yet, we also recognize that there is plenty of room to improve in our model. In particular, more rigorous ways of obtaining information about neural dynamics may be necessary. For instance, our method to estimate stabilized states of a neuronal ensemble may be suboptimal to many different real data. Also, we need to apply our decoder to ensemble data from many different brain regions to generalize it further. Finally, continuous efforts to reduce computational loads in the proposed decoder will be required.


## 4. Conclusions

A number of methods have been continuously developed to decode neuronal ensemble activity. The design of decoding models draws upon the properties of a target brain region, a recording technique, neural signals, and the objectives of decoding. Among them, the maximum entropy decoding method takes into account correlations between neurons as well as firing rates of individual neurons. In particular, it produces good decoding performance when there exist marginal but clear interactions between neurons. However, the current maximum entropy decoder does not capture time-varying characteristics of neuronal ensemble activity, which often deliver essential information about underlying brain functions. Hence, the present study addresses this issue by developing a novel decoder that incorporates dynamical properties of a neuronal ensemble in the model while maintains the key functions of the maximum entropy decoder. We demonstrate that more information can be successfully decoded using the proposed decoder compared to the conventional maximum entropy decoder.



---

*Source: 218373-2014-07-09.xml*