# Structural Credit Assignment with Coordinated Exploration

Stephen Chung University of Massachusetts Amherst, Department of Computer
Science<EMAIL_ADDRESS>
(April 2022)

##  1 Overview

A biologically plausible method for training an Artificial Neural Network
(ANN) involves treating each unit as a stochastic Reinforcement Learning (RL)
agent, thereby considering the network as a team of agents. Consequently, all
units can learn via REINFORCE, a local learning rule modulated by a global
reward signal, which aligns more closely with biologically observed forms of
synaptic plasticity [1, Chapter 15]. However, this learning method tends to be
slow and does not scale well with the size of the network. This inefficiency
arises from two factors impeding effective structural credit assignment: (i)
all units independently explore the network, and (ii) a single reward is used
to evaluate the actions of all units.

Accordingly, methods aimed at improving structural credit assignment can
generally be classified into two categories. The first category includes
algorithms that enable coordinated exploration among units, such as MAP
propagation [2]. The second category encompasses algorithms that compute a
more specific reward signal for each unit within the network, like Weight
Maximization and its variants [3]. In this research report, our focus is on
the first category. We propose the use of Boltzmann machines [4] or a
recurrent network for coordinated exploration. We show that the negative
phase, which is typically necessary to train Boltzmann machines, can be
removed. The resulting learning rules are similar to the reward-modulated
Hebbian learning rule [5]. Experimental results demonstrate that coordinated
exploration significantly exceeds independent exploration in training speed
for multiple stochastic and discrete units based on REINFORCE, even surpassing
straight-through estimator (STE) backpropagation [6].

###  1.1 An Example: How should we vote?

Before we begin, let us consider an intuitive example of how coordinated
exploration may help structural credit assignment. Suppose that you are in a
town of 99 people. At the end of each year, all people in the town vote for a
mayor next year and there are always two candidates in the election: John and
Peter. The candidate with more votes wins and will influence the economy of
the town next year. This is the first year, and people know nothing about the
two candidates. Therefore, you and the others decide whom to vote by flipping
coins. Suppose you voted for John, but Peter won, and Peter performed poorly
as a mayor next year. The question is—whom should you vote for next year?

REINFORCE says that you should increase doing the action that leads to a
positive reward and vice versa. In this case, you should vote more for Peter
the next year because you got a negative reward (the economy was poor next
year) after voting for John. Though it may not seems sensible, you can still
learn the optimal candidate to vote for with this learning method in the long
run—when there are 49 people voting for each person (this happens with the
probability of $\binom{98}{49}0.5^{49}\approx 0.0789$, which decays to $0$
quickly when there are more people in the town) and you can determine the
voting result, you can learn correctly with REINFORCE; other than that, you
are only learning from noise as you have no influences over the election. This
illustrates why REINFORCE does not scale well with the size of networks.

MAP propagation is similar to REINFORCE, but you ‘forgot’ the person you voted
for and ‘believed’ that you have voted for the winning candidate. In this
case, all people in the town, including you, believed that they had voted for
Peter instead. The negative reward is therefore associated with the action of
voting Peter, and all people learn to vote more for John the next year. This
is like coordination in hindsight—all people believed that they had voted for
the same candidate. However, there was actually no coordination between people
when deciding whom to vote, since all people decided whom to vote by flipping
their own coins.

Another form of coordinated exploration is through letting people in the town
communicate with one another before voting. After flipping the coin, you were
prepared to vote for John. But you observed that the majority was prepared to
vote for Peter, and so you followed the majority by voting for Peter. The
negative reward is therefore associated with the action of voting Peter, and
you learn to vote more for John the next year.

This example suggests that coordinated exploration, whether in hindsight or
not, may help with structural credit assignment. We will focus on the second
case, i.e. coordinated exploration without hindsight, in this report.

##  2 Forms of Exploration

For simplicity, we consider an MDP with only immediate reward111The algorithms
in this paper can be applied in general MDPs by replacing the reward with the
sum of discounted reward (i.e. return) or TD error, and can be applied in
supervised learning tasks by replacing the reward with the negative loss.
defined by a tuple $(\mathcal{S},\mathcal{A},R,d_{0})$, where $\mathcal{S}$ is
the set of states, $\mathcal{A}$ is the set of actions,
$R:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$ is the reward function,
and $d_{0}:\mathcal{S}\rightarrow[0,1]$ is the initial state distribution.
Denoting the state, action, and reward by $S$, $A$, and $R$ respectively,
$\operatorname*{\text{Pr}}(S=s)=d_{0}(s)$ and
$\operatorname*{\mathbb{E}}[R|S=s,A=a]=R(s,a)$. We are interested in learning
the policy $\pi:\mathcal{S}\times\mathcal{A}\rightarrow[0,1]$ such that
selecting actions according to $\operatorname*{\text{Pr}}(A=a|S=s)=\pi(s,a)$
maximizes the expected reward $\operatorname*{\mathbb{E}}[R]$.

###  2.1 Independent Exploration

Let us consider the case where the policy is computed by an ANN with one
hidden layer of $N$ Bernoulli-logistic units: to sample action $A$ conditioned
on $S$, we first sample the vector of activation values of the hidden layer,
denoted as $H$, by:

$\operatorname*{\text{Pr}}(H=h|S=s)\propto\exp((W^{T}s+b)^{T}h),$ (1)

where $h\in\\{0,1\\}^{N}$, $W$ is the weight matrix of the hidden layer, and
$b$ is the bias vector of the hidden layer. Equivalently, all $H_{i}$ (we use
subscript $i$ to denote entry $i$ in a vector, and subscript $ij$ to denote
entry $(i,j)$ in a matrix) are independent with each other conditioned on $S$,
and the conditional distribution is given by:

$\operatorname*{\text{Pr}}(H_{i}=h|S=s)=\begin{cases}1-\sigma\left(\sum_{j}W_{ji}s_{j}+b_{i}\right),&\text{if
$h=0$,}\\\ \sigma\left(\sum_{j}W_{ji}s_{j}+b_{i}\right),&\text{if
$h=1$,}\end{cases}$ (2)

where $h\in\\{0,1\\}$ and $\sigma$ denotes the sigmoid function. After
sampling $H$, we sample $A$ from an appropriate distribution conditioned on
$H$. For example, if the set of actions $\mathcal{A}=\\{0,1\\}$, we can use an
output layer of a single Bernoulli-logistic unit:

$\operatorname*{\text{Pr}}(A=a|H=h)=\begin{cases}1-\sigma\left(\sum_{j}W^{out}_{j}h_{j}+b^{out}\right),&\text{if
$a=0$,}\\\ \sigma\left(\sum_{j}W^{out}_{j}h_{j}+b^{out}\right),&\text{if
$a=1$,}\end{cases}$ (3)

where $a\in\\{0,1\\}$, $W^{out}$ is the weight vector of the output layer, and
$b^{out}$ is the scalar bias of the output layer. The parameters to be learned
for this network are: $W,b,W^{out},b^{out}$.

Learning Rule – To derive a learning rule for $W$, we first compute the
gradient of the expected reward $\operatorname*{\mathbb{E}}[R]$ w.r.t. $W$:

$\displaystyle\nabla_{W_{ji}}\operatorname*{\mathbb{E}}[R]$
$\displaystyle=\operatorname*{\mathbb{E}}[R\nabla_{W_{ji}}\log\operatorname*{\text{Pr}}(H|S)]$
(4)
$\displaystyle=\operatorname*{\mathbb{E}}[R(H_{i}-\sigma(\sum_{j}W_{ji}S_{j}+b_{i}))S_{j}],$
(5)
$\displaystyle=\operatorname*{\mathbb{E}}[R(H_{i}-\operatorname*{\mathbb{E}}[H_{i}|S])S_{j}].$
(6)

Therefore, to perform gradient ascent on the expected reward, we can use the
following REINFORCE learning rule [7]:

$\displaystyle W_{ji}$ $\displaystyle\leftarrow W_{ji}+\alpha
R(H_{i}-\operatorname*{\mathbb{E}}[H_{i}|S])S_{j},$ (7)

where $\alpha>0$ denotes the step size. Also note that
$\operatorname*{\mathbb{E}}[H_{i}|S]=\sigma\left(\sum_{j}W_{ji}S_{j}+b_{i}\right)$.

The learning rule is intuitive: the expectation of
$R(H_{i}-\operatorname*{\mathbb{E}}[H_{i}|S])$ is simply the covariance
between $R$ and $H_{i}$ (conditioned on $S$); that is,
$\operatorname*{\mathbb{E}}[R(H_{i}-\operatorname*{\mathbb{E}}[H_{i}|S])|S]=Cov(R,H_{i}|S)$.
Hence if the covariance is positive, the learning rule makes unit $i$ to fire
more, and vice versa.

Recall in statistics that the covariance between two random variables $X$ and
$Y$ equals:

$Cov(X,Y)=\operatorname*{\mathbb{E}}[(X-\operatorname*{\mathbb{E}}[X])(Y-\operatorname*{\mathbb{E}}[Y])]=\operatorname*{\mathbb{E}}[X(Y-\operatorname*{\mathbb{E}}[Y])]=\operatorname*{\mathbb{E}}[(X-\operatorname*{\mathbb{E}}[X])Y],$
(8)

and thus the following three estimators that estimate the covariance are all
unbiased: (i)
$(X-\operatorname*{\mathbb{E}}[X])(Y-\operatorname*{\mathbb{E}}[Y])$, (ii)
$X(Y-\operatorname*{\mathbb{E}}[Y])$, (iii)
$(X-\operatorname*{\mathbb{E}}[X])Y$. We can further show that the first
estimator has the lowest variance among the three. In other words, an
estimator with two-sided centering gives a lower variance than an estimator
with one-sided centering, and an estimator with one-sided centering is already
an unbiased estimator of the covariance.

Noting that (7) only uses one-sided centering to estimate the covariance, an
alternative learning rule with the same expected update is to use two-sided
centering:

$\displaystyle W_{ji}$ $\displaystyle\leftarrow
W_{ji}+\alpha(R-\operatorname*{\mathbb{E}}[R|S])(H_{i}-\operatorname*{\mathbb{E}}[H_{i}|S])S_{j},$
(9)

which is equivalent to REINFORCE with $\operatorname*{\mathbb{E}}[R|S]$ as
baseline [7]. $\operatorname*{\mathbb{E}}[R|S]$ is generally unknown and so
needs to be estimated by another network which is usually called a critic
network.

Finally, we can also do one-sided centering on $R$ only, leading to the
following learning rule that has the same expected update:

$\displaystyle W_{ji}$ $\displaystyle\leftarrow
W_{ji}+\alpha(R-\operatorname*{\mathbb{E}}[R|S])H_{i}S_{j},$ (10)

which has the property that the parameter of a unit is not updated if the unit
does not fire—when $H_{i}=0$, the above updates equal zero. This property may
be desirable as biological learning rules such as reward-modulated spike-
timing-dependent plasticity (R-STDP) [5] also require a neuron to fire to
affect synaptic strength.

All three learning rules (7), (9) and (10) follow gradient ascent on
$\operatorname*{\mathbb{E}}[R]$ in expectation. Among the three learning
rules, learning rule (9) is associated with the lowest variance due to the use
of two-sided centering. However, as we will see in the next section, learning
rule (10) can be generalized more easily to incorporate coordinated
exploration.

Note that a similar discussion applies to the bias $b$, and the learning rules
of $b_{i}$ are the same as above but with $S_{j}$ set to $1$.

###  2.2 Coordinated Exploration with Boltzmann Machines

As discussed previously, independent exploration coupled with a uniform reward
signal makes learning inefficient and scales poorly with the number of units
in the network. In this section we consider how to generalize the above
network to have coordinated exploration.

Let us consider the conditional distribution of $H$ given by (1), which says
that hidden units are conditionally independent with one another. A natural
way to allow interaction between hidden units is to add a cross term between
hidden units:

$\operatorname*{\text{Pr}}(H=h|S=s)\propto\exp((W^{T}s+b)^{T}h+ch^{T}W^{rec}h),$
(11)

where $c\geq 0$ is a scalar hyperparameter that controls the strength of
interaction between hidden units, and $W^{rec}$ is a symmetric weight matrix
with a zero diagonal that determines the interaction dynamics between hidden
units. We call $W^{rec}$ the recurrent weight and $W$ the feedforward weight
to distinguish between these two weights. Note that we recover the case of
independent exploration in (1) by setting $c=0$. The distribution given by
(11) can be viewed as a fully connected Boltzmann machine [4], with $W^{rec}$
being the weight and $W^{T}s+b$ being the bias of the Boltzmann machine.

Figure 1: Illustration of how $H$ is generated by independent exploration,
coordinated exploration with Boltzmann machines and recurrent networks. Left
figure: a standard one-hidden-layer network of Bernoulli-logistic units. Each
$H_{i}$ is sampled from $Ber((W^{T}S+b)_{i})$. Middle figure: the distribution
of $H$ is given by a Boltzmann machine with weight $W^{rec}$ and bias
$W^{T}S+b$. Right figure: the hidden layer is a recurrent layer with recurrent
weight $W^{rec}$ and common input $W^{T}S+b$, so each $H_{(t+1),i}$ is sampled
from $Ber((c(W^{rec})^{T}H_{(t)}+W^{T}S+b)_{i})$. If $W^{rec}$ is symmetric
and has a zero diagonal, as $T\rightarrow\infty$, the distribution of $H$
depicted in the middle figure and that of $H_{(T)}$ depicted in the right
figure become the same (to be rigorous, only one random entry of $H_{(t)}$ is
updated in each recurrent step in the right figure for this to hold as
asynchronous sampling is required).

We can then compute the distribution of $H_{i}$ conditioned on $H_{-i}$ and S
($H_{-i}$ denotes all entries of $H$ except $H_{i}$):

$\operatorname*{\text{Pr}}(H_{i}=h_{i}|H_{-i}=h_{-i},S=s)=\begin{cases}1-\sigma\left(c\sum_{j}W^{rec}_{ji}h_{j}+\sum_{j}W_{ji}s_{j}+b_{i}\right),&\text{if
$h_{i}=0$,}\\\
\sigma\left(c\sum_{j}W^{rec}_{ji}h_{j}+\sum_{j}W_{ji}s_{j}+b_{i}\right),&\text{if
$h_{i}=1$.}\end{cases}$ (12)

To sample $H$, we can perform Gibbs sampling: we first initialize $H_{(0)}$ to
some values such as a zero vector (we use the subscript $(t)$ to denote time
step $t$ as opposed to entry), then we randomly select one unit and replace
its activation value by sampling from (12) to obtain $H_{(1)}$. We continue
this way to obtain $H_{(2)},H_{(3)},...$, and it can be shown that the
limiting distribution of $H_{(t)}$ as $t\rightarrow\infty$ is the same as the
distribution given by (11).

However, this asynchronous sampling method is computationally
expensive222However, asynchronous computation should be closer to the dynamics
of biological neurons as biological neurons do not fire at fixed intervals..
Therefore, we use synchronous sampling method instead, that is, we update the
activation values of all hidden units on the same time step:

$H_{(t+1),i}|H_{(t)},S\sim
Ber\left(\sigma\left(c\sum_{j}W^{rec}_{ji}H_{(t),j}+\sum_{j}W_{ji}S_{j}+b_{i}\right)\right).$
(13)

This allows easy sampling of $H_{(t+1)}$ conditioned on $H_{(t)}$ and $S$.
Then we let $H=H_{(T)}$, where $T\geq 1$ is a hyperparameter determining the
number of steps for Gibbs sampling, and $H$ is passed to the next layer as
input. See Figure 1 for an illustration.

Though the limiting distribution of $H_{(t)}$ as $t\rightarrow\infty$ is not
guaranteed to converge due to synchronous sampling, we found that it works
well in experiments. We also found in experiments that $T\geq 25$ works quite
well.

Learning Rule – We now consider how to train a network when the distribution
of $H$ is given by (11). Let us first compute
$\nabla_{W^{rec}}\log\operatorname*{\text{Pr}}(H=h|S=s)$, which will be useful
when deriving learning rules later on:

$\displaystyle\nabla_{W^{rec}_{ji}}\log\operatorname*{\text{Pr}}(H=h|S=s)$
(14) $\displaystyle=$
$\displaystyle\nabla_{W^{rec}_{ji}}\log\frac{\exp((W^{T}s+b)^{T}h+ch^{T}W^{rec}h)}{\sum_{h}\exp((W^{T}s+b)^{T}h+ch^{T}W^{rec}h)}$
(15) $\displaystyle=$ $\displaystyle
ch_{i}h_{j}-\nabla_{W^{rec}_{ji}}\log\sum_{h}\exp((W^{T}s+b)^{T}h+ch^{T}W^{rec}h)$
(16) $\displaystyle=$ $\displaystyle
ch_{i}h_{j}-\frac{\sum_{h}\exp((W^{T}s+b)^{T}h+ch^{T}W^{rec}h)\nabla_{W^{rec}_{ji}}((W^{T}s+b)^{T}h+ch^{T}W^{rec}h)}{\sum_{h}\exp((W^{T}s+b)^{T}h+ch^{T}W^{rec}h)}$
(17) $\displaystyle=$ $\displaystyle
ch_{i}h_{j}-c\sum_{h}\operatorname*{\text{Pr}}(H=h|S=s)h_{i}h_{j}$ (18)
$\displaystyle=$ $\displaystyle
c(h_{i}h_{j}-\operatorname*{\mathbb{E}}[H_{i}H_{j}|S=s]).$ (19)

This derivation is similar to the one used to derive the learning rule of
Contrastive Divergence for Boltzmann machines. The term
$\operatorname*{\mathbb{E}}[H_{i}H_{j}|S=s]$ is called negative statistics and
is usually estimated by drawing samples from the Boltzmann machine in the
unclamped phase.

Similarly, we can show that:

$\displaystyle\nabla_{W_{ji}}\log\operatorname*{\text{Pr}}(H=h|S=s)$
$\displaystyle=c(h_{i}-\operatorname*{\mathbb{E}}[H_{i}|S=s])s_{j},$ (20)
$\displaystyle\nabla_{b_{i}}\log\operatorname*{\text{Pr}}(H=h|S=s)$
$\displaystyle=c(h_{i}-\operatorname*{\mathbb{E}}[H_{i}|S=s]).$ (21)

Now we are ready to compute the gradient of the expected reward w.r.t. the
parameters. We first consider the gradient of the expected reward w.r.t.
$W^{rec}$:

$\displaystyle\nabla_{W^{rec}_{ji}}\operatorname*{\mathbb{E}}[R]$
$\displaystyle=\operatorname*{\mathbb{E}}[R\nabla_{W^{rec}_{ji}}\log\operatorname*{\text{Pr}}(H|S)]$
(22)
$\displaystyle=\operatorname*{\mathbb{E}}[cR(H_{i}H_{j}-\operatorname*{\mathbb{E}}[H_{i}H_{j}|S])].$
(23)

The last line follows from (19). The gradient w.r.t other parameters are
similar:

$\displaystyle\nabla_{W_{ji}}\operatorname*{\mathbb{E}}[R]$
$\displaystyle=\operatorname*{\mathbb{E}}[cR(H_{i}-\operatorname*{\mathbb{E}}[H_{i}|S])S_{j}],$
(24) $\displaystyle\nabla_{b_{i}}\operatorname*{\mathbb{E}}[R]$
$\displaystyle=\operatorname*{\mathbb{E}}[cR(H_{i}-\operatorname*{\mathbb{E}}[H_{i}|S])].$
(25)

Therefore, to perform gradient ascent on the expected reward, we can use the
following learning rules:

$\displaystyle W^{rec}_{ji}$ $\displaystyle\leftarrow W^{rec}_{ji}+\alpha
cR(H_{i}H_{j}-\operatorname*{\mathbb{E}}[H_{i}H_{j}|S]),$ (26) $\displaystyle
W_{ji}$ $\displaystyle\leftarrow W_{ji}+\alpha
R(H_{i}-\operatorname*{\mathbb{E}}[H_{i}|S])S_{j},$ (27) $\displaystyle b_{i}$
$\displaystyle\leftarrow b_{i}+\alpha
R(H_{i}-\operatorname*{\mathbb{E}}[H_{i}|S]).$ (28)

Note that (27) is the same as (7), but the expected value
$\operatorname*{\mathbb{E}}[H_{i}|S]$ has to be computed or estimated
differently since the distributions of $H$ are not the same.

These learning rules require the computation or estimation of the negative
statistics $\operatorname*{\mathbb{E}}[H_{i}H_{j}|S]$ and
$\operatorname*{\mathbb{E}}[H_{i}|S]$. Both terms are intractable to compute
directly, and so we have to resort to Monte Carlo estimation as in Boltzmann
machines: we can perform multiple Gibbs samplings to obtain samples of $H$
conditioned on $S$, and average these samples to estimate the negative
statistics. But this is hardly biologically plausible and also computationally
inefficient. Is it possible to remove the negative statistics from the
learning rules while still performing gradient ascent?

Let us consider the terms in the learning rules.
$R(H_{i}H_{j}-\operatorname*{\mathbb{E}}[H_{i}H_{j}|S])$ is estimating the
covariance between $R$ and $H_{i}H_{j}$, and
$R(H_{i}-\operatorname*{\mathbb{E}}[H_{i}|S])$ is estimating the covariance
between $R$ and $H_{i}$. Thus, the discussion in the previous section on
different estimators for covariance can also be applied here. In particular,
we can do one-sided centering on $R$ instead of $H$ to obtain the following
learning rules with the same expected update:

$\displaystyle W^{rec}_{ji}$ $\displaystyle\leftarrow W^{rec}_{ji}+\alpha
c(R-\operatorname*{\mathbb{E}}[R|S])H_{i}H_{j},$ (29) $\displaystyle W_{ji}$
$\displaystyle\leftarrow
W_{ji}+\alpha(R-\operatorname*{\mathbb{E}}[R|S])H_{i}S_{j},$ (30)
$\displaystyle b_{i}$ $\displaystyle\leftarrow
b_{i}+\alpha(R-\operatorname*{\mathbb{E}}[R|S])H_{i}.$ (31)

These learning rules follow gradient ascent on $\operatorname*{\mathbb{E}}[R]$
in expectation. Again, we require a critic network to estimate
$\operatorname*{\mathbb{E}}[R|S]$, but this is considerably simpler than
estimating the negative statistics. In fact,
$R-\operatorname*{\mathbb{E}}[R|S]$ is simply TD error in a single-time-step
MDP, and can be replaced with TD error in a multiple-time-step MDP, noting
that TD error is also centered. In this way, we can get rid of the negative
statistics in the learning rule.

It should be noted that we need at least one-side centering for learning to
occur. Consider the learning rule with no centering: $W^{rec}_{ji}\leftarrow
W^{rec}_{ji}+\alpha cRH_{i}H_{j}$; if the reward signals from the environment
are all non-negative (e.g. Cartpole), then the weight can only be increased
but not decreased, and learning cannot occur.

The learning rules are similar to reward-modulated Hebbian learning [5]—for
example, (29) is the product between three terms: (i) TD error, (ii) the pre-
synaptic activation value $H_{j}$ and (iii) the post-synaptic activation value
$H_{i}$. Furthermore, as long as $W^{rec}$ is initialized to be symmetric, it
will stay symmetric during training since the updates to $W^{rec}_{ji}$ and
$W^{rec}_{ij}$ are the same, so we do not need additional mechanisms to ensure
symmetric connections as in backprop.

1Algorithm Hyperparameter: Number of sampling steps $T\geq 1$; strength of
recurrent connections $c\geq 0$; number of hidden units $N\geq 1$; step size
$\alpha>0$;

2 Initialize parameters: $W\in\mathbb{R}^{dim(S)\times
N},b\in\mathbb{R}^{N},W^{rec}\in\mathbb{R}^{N\times N}$ (symmetric with zero
diagonal), parameter for the output layer $\theta^{out}$ ;

3 Loop forever:  __

4 Sample state $S$ ;

/* Sampling $H$ and $A$ */

5 Sample $H_{i}\sim
Ber\left(\sigma\left(\sum_{j}W_{ji}S_{j}+b_{i}\right)\right)$ for all
$i\in\\{1,2,...,N\\}$;

6 for  _$t:=1,2,...,T$_ do

7 Sample $H^{\prime}_{i}\sim
Ber\left(\sigma\left(c\sum_{j}W^{rec}_{ji}H_{j}+\sum_{j}W_{ji}S_{j}+b_{i}\right)\right)$
for all $i\in\\{1,2,...,N\\}$ ;

8 $H\leftarrow H^{\prime}$ ;

9

10 Sample $A\sim\operatorname*{\text{Pr}}(A=\cdot|H;\theta^{out})$ ;

11 Take action $A$, observe $R$ and receive estimate
$\hat{\operatorname*{\mathbb{E}}}[R|S]$ from critic network;

/* REINFORCE */

12 $W^{rec}_{ji}\leftarrow W^{rec}_{ji}+\alpha
c(R-\hat{\operatorname*{\mathbb{E}}}[R|S])H_{i}H_{j}$ for all $i\neq
j,i,j\in\\{1,2,...,N\\}$ ;

13 $W_{ji}\leftarrow
W_{ji}+\alpha(R-\hat{\operatorname*{\mathbb{E}}}[R|S])H_{i}S_{j}$ for all
$i\in\\{1,2,...,N\\},j\in\\{1,2,...,dim(S)\\}$ ;

14 $b_{i}\leftarrow
b_{i}+\alpha(R-\hat{\operatorname*{\mathbb{E}}}[R|S])H_{i}$ for all
$i\in\\{1,2,...,N\\}$;

15
$\theta^{out}\leftarrow\theta^{out}+\alpha(R-\hat{\operatorname*{\mathbb{E}}}[R|S])\nabla_{\theta^{out}}\log
P(A|H)$;

16 Train the critic network with error signal
$R-\hat{\operatorname*{\mathbb{E}}}[R|S]$;

17

Algorithm 1 Coordinated Exploration with a Boltzmann Machine Figure 2:
Illustration of learning by coordinated exploration with Boltzmann machines.
First, we initialize $H_{(0)}$ and then minimize the energy, $E(H)=-\log
P(H|S)$, by updating $H$ stochastically to obtain $H_{(T)}$, as shown in the
left figure. Then, we pass $H_{(T)}$ to the next layer and receive the TD
error. Finally, if the TD error is positive, we decrease the energy at
$H_{(T)}$ by updating the parameters of the network, as shown by the shift
from the black curve to the blue dotted curve in the middle figure. This makes
$H$ easier to ‘fall’ into the same location in future sampling. The case for
negative TD error is the opposite and is shown in the right figure.

The full algorithm is given in Algorithm 1. An illustration of the algorithm
from the perspective of energy-based learning is shown in Figure 2. Note that
if we set $c$, the strength of recurrent connections, to zero, then we recover
the case of independent exploration trained with REINFORCE (one-sided
centering on $R$) given by the learning rule (10).

Experiments – To test the algorithm, we consider the $k$-bit multiplexer task.
The input is sampled from all possible values of a binary vector of size
$k+2^{k}$ with equal probability. The action set is $\\{0,1\\}$, and we give a
reward of $+1$ if the action equals the multiplexer’s output and $-1$
otherwise. We consider $k=4$ here, so the dimension of the state space is 20
and there are $2^{20}=1048576$ possible states. We call a single interaction
of the network from receiving the state to receiving the reward an episode;
the series of episodes as the network’s performance increases during training
is called a run.

The update step size $\alpha$ for gradient ascent is chosen as $0.005$, and
the batch size is chosen as $16$ (i.e. we estimated the gradient for $16$
episodes and averaged them in the update). We used the Adam optimizer [8].
These hyperparameters are chosen prior to the experiments without any tuning.
We use a one-hidden-layer ANN as the critic network and train it by backprop.
These settings are the same for all our experiments.

First, we tested Algorithm 1 with $N=64$ hidden units and $T=25$ sampling
steps. The learning curve (i.e. the reward of each episode throughout a run)
with varying $c$ (the strength of recurrent connection) is shown in Figure 3.
The figure also shows the learning curve of a network with one hidden layer of
$64$ Bernoulli-logistic units trained by REINFORCE with baseline (9), which we
label ‘REINFORCE (indep)’. We observe that:

  * •

Both ‘REINFORCE (indep)’ and ‘Alg. 1; c=0.00’ use independent exploration. The
former one (9) use centering on both side while the latter one (10) only use
centering on $R$. As expected, two-sided centering gives a lower variance and
so the former one performs better.

  * •

However, as $c$ increases from $0$ to $0.25$, the performance of Algorithm 1
improves due to the coordinated exploration, offsetting the reduced
performance from one-sided centering.

  * •

The optimal $c$ is $0.25$; at $c=0.25$, the learning curve, especially during
the first $1e6$ episodes, is significantly better than ‘REINFORCE (indep)’,
indicating that Algorithm 1 can improve learning speed significantly with a
proper $c$.

  * •

But as $c$ increases larger than $0.25$, performance deteriorates since inputs
from other hidden units instead of the state dominate the unit. In particular,
the asymptotic performance is unstable for $c\geq 0.75$.

The experiment result suggests that coordinated exploration with Boltzmann
machines can improve learning speed; the benefit from coordinated exploration
offsets the reduced performance from one-sided centering. Also, $c$, the
strength of recurrent connection, is an important hyperparameter that needs to
be carefully tuned. Simply setting $c=1$ can lead to reduced performance.

Figure 3: Learning curve of Algorithm 1 with varying $c$ and $T=25$. Results
are averaged over 5 independent runs, and the error bar represents the
standard deviation over the runs.

To understand how the different learning rules and networks scale with $N$,
the number of hidden units, we repeat the above experiments for
$N\in\\{8,16,32,64,96,128\\}$. We used $T=25$, and $c$ was tuned for each $N$
to maximize the average reward throughout all episodes. The average reward
during the first 1e6 episodes (to evaluate the learning speed) and the last
1e6 episodes (to evaluate the asymptotic performance) for different $N$ is
shown in Figure 4 and 5 respectively. For comparison, we also show the results
from training a network with one hidden layer of $N$ Bernoulli-logistic units
using REINFORCE with baseline (9), unbiased Weight Maximization, and STE
backprop.

From Figure 4, we observe that Algorithm 1 has the best learning speed among
all algorithms considered, and also scales well with the number of hidden
units. It is worth noting that _Algorithm 1 is faster than STE backprop_ at
all $N$ tested. From Figure 5, we observe that the asymptotic performance of
Algorithm 1 is similar to other algorithms considered except Algorithm 2,
another algorithm that will be discussed in the next section. However, it
should be noted that the computation cost of Algorithm 1 is more expensive
than other algorithms considered. As we use $T=25$, it means that we require
$25$ times computational cost to sample $H$ compared to other algorithms.

Figure 4: Average reward during the first 1e6 episodes for a network with
varying $N$. Results are averaged over 5 independent runs, and the error bar
represents the standard deviation over the runs. Figure 5: Average reward
during the last 1e6 episodes for a network with varying $N$. Results are
averaged over 5 independent runs, and the error bar represents the standard
deviation over the runs.

Additional experiments, such as ablation analysis on $T$ and how the optimal
$c$ changes with $N$, can be found in Appendix A.1.

###  2.3 Coordinated Exploration with Recurrent Networks

Though symmetric recurrent weights arise naturally from Algorithm 1, symmetric
connections between biological neurons are not commonly observed. Learning
rules such as R-STDP is asymmetric instead of symmetric: if the post-synaptic
spike arrives before the pre-synaptic spike arrives, then the synaptic
strength is decreased instead of increased. Moreover, the use of one-sided
centering increases variance and hurts performance. Is it possible to remove
the requirement of symmetric connection in Algorithm 1 while using two-sided
centering?

Let look at how $H$ is generated in Algorithm 1, as given by (13). One can
also see this equation as a recurrent neural network, with $S$ being the
common input at all time steps, as shown in Figure 1. From this perspective,
we can simply apply REINFORCE on each time step to obtain an unbiased
estimator of the gradient. For example, the gradient of the expected reward
w.r.t. $W^{rec}$ can be computed as:

$\displaystyle\nabla_{W^{rec}_{ji}}\operatorname*{\mathbb{E}}[R]$
$\displaystyle=\operatorname*{\mathbb{E}}\left[R\sum_{t=1}^{T}\nabla_{W^{rec}_{ji}}\log\operatorname*{\text{Pr}}(H_{(t+1)}|S,H_{(t)})\right]$
(32)
$\displaystyle=\operatorname*{\mathbb{E}}\left[cR\sum_{t=1}^{T}(H_{(t+1),i}-\operatorname*{\mathbb{E}}[H_{(t+1),i}|S,H_{(t)}])H_{(t),j}\right],$
(33)

where
$\operatorname*{\mathbb{E}}[H_{(t+1),i}|S,H_{(t)}]=\sigma\left(c\sum_{j}W^{rec}_{ji}H_{(t),j}+\sum_{j}W_{ji}S_{j}+b_{i}\right)$.
The first line follows by the fact that a recurrent ANN can be viewed as a
multi-layer ANN with shared weights, and so we can apply REINFORCE on each
layer.

The gradient w.r.t. $W$ and $b$ can be computed similarly:

$\displaystyle\nabla_{W_{ji}}\operatorname*{\mathbb{E}}[R]$
$\displaystyle=\operatorname*{\mathbb{E}}\left[R\sum_{t=1}^{T}(H_{(t+1),i}-\operatorname*{\mathbb{E}}[H_{(t+1),i}|S,H_{(t)}])S_{j}\right],$
(34) $\displaystyle\nabla_{b_{i}}\operatorname*{\mathbb{E}}[R]$
$\displaystyle=\operatorname*{\mathbb{E}}\left[R\sum_{t=1}^{T}(H_{(t+1),i}-\operatorname*{\mathbb{E}}[H_{(t+1),i}|S,H_{(t)}])\right].$
(35)

Note that these equations do not assume $W^{rec}$ to be symmetric or to have
zero diagonal. It also does not require $H^{(T)}$ to converge to any limiting
distribution, so $T$ can be small. From these equations, we arrive at the
following learning rules that follow gradient ascent on expected reward in
expectation:

$\displaystyle W^{rec}_{ji}$ $\displaystyle\leftarrow W^{rec}_{ji}+\alpha
cR\sum_{t=1}^{T}(H_{(t+1),i}-\operatorname*{\mathbb{E}}[H_{(t+1),i}|S,H_{(t)}])H_{(t),j},$
(36) $\displaystyle W_{ji}$ $\displaystyle\leftarrow W_{ji}+\alpha
R\sum_{t=1}^{T}(H_{(t+1),i}-\operatorname*{\mathbb{E}}[H_{(t+1),i}|S,H_{(t)}])S_{j},$
(37) $\displaystyle b_{i}$ $\displaystyle\leftarrow b_{i}+\alpha
R\sum_{t=1}^{T}(H_{(t+1),i}-\operatorname*{\mathbb{E}}[H_{(t+1),i}|S,H_{(t)}]).$
(38)

However, these learning rules suffer from a large variance as they can be
viewed as REINFORCE applied to a multi-layer ANN with shared weights. To solve
this issue, we borrow the idea of temporal credit assignment in RL. Consider
that we use a very large $T$, e.g. $T=1000$. It is unlikely that $H_{(t)}$ for
a small $t$ can have influences on $H_{(1000)}$ and also the reward $R$. In
other words, the dependency between $R$ and $H_{(t)}$ decreases as $t$
decreases from $T$. Thus, instead of placing equal importance on every time
step in the learning rule, we can multiply the terms in the learning rule with
an exponentially decaying rate $0\leq\lambda\leq 1$:

$\displaystyle W^{rec}_{ji}$ $\displaystyle\leftarrow W^{rec}_{ji}+\alpha
cR\sum_{t=1}^{T}\lambda^{T-t}(H_{(t+1),i}-\operatorname*{\mathbb{E}}[H_{(t+1),i}|S,H_{(t)}])H_{(t),j},$
(39) $\displaystyle W_{ji}$ $\displaystyle\leftarrow W_{ji}+\alpha
R\sum_{t=1}^{T}\lambda^{T-t}(H_{(t+1),i}-\operatorname*{\mathbb{E}}[H_{(t+1),i}|S,H_{(t)}])S_{j},$
(40) $\displaystyle b_{i}$ $\displaystyle\leftarrow b_{i}+\alpha
R\sum_{t=1}^{T}\lambda^{T-t}(H_{(t+1),i}-\operatorname*{\mathbb{E}}[H_{(t+1),i}|S,H_{(t)}]).$
(41)

These learning rules can be easily implemented online by eligibility traces.
For example, the learning rule for $W^{rec}$ is equivalent to:

$\displaystyle\mathbf{z}^{rec}_{(0),ji}$ $\displaystyle=0.$ (42)
$\displaystyle\mathbf{z}^{rec}_{(t+1),ji}$
$\displaystyle=\lambda\mathbf{z}^{rec}_{(t),ji}+(H_{(t+1),i}-\operatorname*{\mathbb{E}}[H_{(t+1),i}|S,H_{(t)}])H_{(t),j},\text{for
}t\in\\{1,2,...,T\\},$ (43) $\displaystyle W^{rec}_{ji}$
$\displaystyle\leftarrow\alpha cR\mathbf{z}^{rec}_{(T),ji}.$ (44)

1Algorithm Hyperparameter: Number of sampling steps $T\geq 1$; strength of
recurrent connections $c\geq 0$; number of hidden units $N\geq 1$; eligibility
decay rate $0\leq\lambda\leq 1$; step size $\alpha>0$;

2 Initialize parameters: $W\in\mathbb{R}^{dim(S)\times
N},b\in\mathbb{R}^{N},W^{rec}\in\mathbb{R}^{N\times N}$, parameter for the
output layer $\theta^{out}$ ;

3 Loop forever:  __

4 Sample state $S$ ;

/* Sampling $H$ and $A$ */

5 Sample $H_{i}\sim
Ber\left(\sigma\left(\sum_{j}W_{ji}S_{j}+b_{i}\right)\right)$ for all
$i\in\\{1,2,...,N\\}$;

6 Initialize $\mathbf{z}^{rec}=\mathbf{0},\mathbf{z}=\mathbf{0}$;

7 for  _$t:=1,2,...,T$_ do

8
$p_{i}\leftarrow\sigma\left(c\sum_{j}W^{rec}_{ji}H_{j}+\sum_{j}W_{ji}S_{j}+b_{i}\right)$
for all $i\in\\{1,2,...,N\\}$ ;

9 Sample $H^{\prime}_{i}\sim Ber(p_{i})$ for all $i\in\\{1,2,...,N\\}$ ;

10
$\mathbf{z}^{rec}_{ji}\leftarrow\lambda\mathbf{z}^{rec}_{ji}+(H^{\prime}_{i}-p_{i})H_{j}$
for all $i,j\in\\{1,2,...,N\\}$ ;

11 $\mathbf{z}_{i}\leftarrow\lambda\mathbf{z}_{i}+(H^{\prime}_{i}-p_{i})$ for
all $i\in\\{1,2,...,N\\}$;

12 $H\leftarrow H^{\prime}$ ;

13

14 Sample $A\sim\operatorname*{\text{Pr}}(A=\cdot|H;\theta^{out})$ ;

15 Take action $A$, observe $R$ and receive estimate
$\hat{\operatorname*{\mathbb{E}}}[R|S]$ from critic network;

/* REINFORCE */

16 $W^{rec}_{ji}\leftarrow W^{rec}_{ji}+\alpha
c(R-\hat{\operatorname*{\mathbb{E}}}[R|S])\mathbf{z}^{rec}_{ji}$ for all
$i,j\in\\{1,2,...,N\\}$ ;

17 $W_{ji}\leftarrow
W_{ji}+\alpha(R-\hat{\operatorname*{\mathbb{E}}}[R|S])\mathbf{z}_{i}S_{j}$ for
all $i\in\\{1,2,...,N\\},j\in\\{1,2,...,dim(S)\\}$ ;

18 $b_{i}\leftarrow
b_{i}+\alpha(R-\hat{\operatorname*{\mathbb{E}}}[R|S])\mathbf{z}_{i}$ for all
$i\in\\{1,2,...,N\\}$;

19
$\theta^{out}\leftarrow\theta^{out}+\alpha(R-\hat{\operatorname*{\mathbb{E}}}[R|S])\nabla_{\theta^{out}}\log
P(A|H)$;

20 Train the critic network with error signal
$R-\hat{\operatorname*{\mathbb{E}}}[R|S]$;

21

Algorithm 2 Coordinated Exploration with a Recurrent Network

The full algorithm is given in Algorithm 2. We also adjust the reward by the
expected reward from a critic network in the algorithm (though this is not
necessary for learning to occur, as opposed to Algorithm 1). It should be
noted that at $c=0$ and $\lambda=0$, the algorithm recovers REINFORCE with
baseline (9) applied on a network with one hidden layer of $N$ Bernoulli-
logistic units.

In addition, at $\lambda=0$ (which works well in our experiments), that is,
the hidden units only learn at the last time step, Algorithm 2 becomes the
same as Algorithm 1, except: (i) Algorithm 2 uses two-sided centering as
$\operatorname*{\mathbb{E}}[H_{(T)}|S,H_{(T-1)}]$ is subtracted from
$H_{(T)}$, which makes the learning rule no longer symmetric and so is the
recurrent weight $W^{rec}$; (ii) the pre-synaptic signal in the learning rule
of the recurrent weight $W^{rec}$ is $H_{(T-1)}$ in Algorithm 2 instead of
$H_{(T)}$ in Algorithm 1. In other words, the pre-synaptic signal is one time
step before the post-synaptic signal in the learning rule of $W^{rec}$ in
Algorithm 2, making it closer to R-STDP instead of reward-modulated Hebbian
learning.

Though Algorithm 2 may seem similar to Algorithm 1 as $H$ is generated in the
exact same way, the dynamics of these two algorithms are very different
throughout learning. Algorithm 2 treats the hidden layer as a recurrent layer
instead of a Boltzmann machine, so any $T\geq 1$ already gives a good
performance instead of $T\geq 25$ in Algorithm 1. The asymmetric connection
also makes $H$ not guaranteed to converge to any distributions as in Boltzmann
machines. These different dynamics manifest in our experiments, which will be
discussed in the following section.

Experiments – To test the algorithm, we use the same setting as in Section
2.2; that is, we consider the task of $4$-bit multiplexer tasks with the same
set of hyperparameters. First, we tested Algorithm 2 with $N=64$ hidden units,
$T=2$ sampling steps and $\lambda=0.25$ trace decay rate. The learning curve
with varying $c$ (the strength of recurrent connection) is shown in Figure 6.
The figure also shows the learning curve of a network with one hidden layer of
$64$ Bernoulli-logistic units trained by REINFORCE with baseline (9), which we
label ‘REINFORCE (indep)’. The learning curve of Algorithm 1 with $c=0.25$
from Figure 3 is also shown here for comparison. We observe that:

Figure 6: Learning curve of Algorithm 2 with varying $c$, $T=2$ and
$\lambda=0.25$. Results are averaged over 5 independent runs, and the error
bar represents the standard deviation over the runs.

  * •

Both ‘REINFORCE (indep)’ and ‘Alg. 2; c=0.00’ use independent exploration, and
the former one is equivalent to Algorithm 2 with both $c=0$ and $\lambda=0$.
Since ‘Alg. 2; c=0.00’ uses $\lambda=0.25$ and there are no interaction
between hidden units, the use of $\lambda$ only adds noise to the parameter
update and leads to a slightly worse performance.

  * •

However, as $c$ increases from $0$ to $0.25$, the performance of Algorithm 2
improves due to the coordinated exploration.

  * •

The optimal $c$ is $0.25$; at $c=0.25$, the learning curve is better than
‘REINFORCE (indep)’ in terms of both learning speed and asymptotic
performance; compared to ‘Alg. 1; c=0.25’, the learning speed is lower, but
the asymptotic performance is better.

  * •

But as $c$ increases larger than $0.25$, performance deteriorates since inputs
from other hidden units instead of the state dominate the unit.

  * •

The asymptotic performance is stable for all $c$, unlike Algorithm 1.

The experiment result suggests that coordinated exploration with recurrent
networks can improve learning speed. Again, $c$, the strength of recurrent
connection, is an important hyperparameter that needs to be carefully tuned.
The dynamics of Algorithm 1 and 2 are also different; Algorithm 1 has a high
learning speed while Algorithm 2 has a better and stable asymptotic
performance. In addition, Algorithm 1 requires $T\geq 25$ for good performance
while the performance of Algorithm 2 remains stable for all $T\geq 1$.

To understand how the different learning rules and networks scale with $N$,
the number of hidden units, we repeat the above experiments for
$N\in\\{8,16,32,64,96,128\\}$. We used $T=2$, and $c$ was tuned for each $N$
to maximize the average reward throughout all episodes. The average reward
during the first 1e6 episodes (to evaluate the learning speed) and the last
1e6 episodes (to evaluate the asymptotic performance) for different $N$ is
shown in Figure 4 and 5 respectively.

From Figure 4, we observe that the learning speed of Algorithm 2 is moderately
higher than ‘REINFORCE (indep)’ and also does not benefit much from increasing
the network’s size. From Figure 5, however, we observe that the asymptotic
performance of Algorithm 2 is the best among all algorithms considered at all
$N$.

The lower learning speed but better asymptotic performance can be explained by
the assumption used in deriving the learning rules. In Algorithm 2, we view
$H$ as the output of a recurrent layer, and we apply REINFORCE to the hidden
units on the last few steps (due to eligibility traces). This is equivalent to
training an ANN of $T$ hidden layers with shared weight $W^{rec}$ by applying
REINFORCE on the last few layers. Thus, the low learning speed of REINFORCE in
training an ANN also manifests here; however, the deeper network increases the
network’s capacity, thereby giving a better asymptotic performance.

Additional experiments, such as ablation analysis on $T$ and $\lambda$, and
how the optimal $c$ changes with $N$, can be found in Appendix A.2.

###  2.4 Discussion

In this report we propose and discuss two different algorithms to allow more
efficient structural credit assignment by coordinated exploration. The first
algorithm generalizes hidden layers of Bernoulli-logistic units to Boltzmann
machines, while the second algorithm generalizes hidden layers to recurrent
layers. It turns out that the activation values of hidden layers in these two
algorithms can be sampled in the exact same way. However, the resulting
learning rules are different and lead to different dynamics. Coordinated
exploration with Boltzmann machines (Algorithm 1) increases learning speed
significantly such that the learning speed is even higher than STE backprop,
whereas coordinated exploration with recurrent networks (Algorithm 2) improves
asymptotic performance. It remains to be seen whether we can combine the
advantages of both algorithms in a single learning rule. Nonetheless, in
either algorithm, we observe an improved learning speed compared to the
baseline of training a network of Bernoulli-logistic units by REINFORCE,
indicating that coordinated exploration can facilitate structural credit
assignment.

In addition, coordinated exploration by Boltzmann machines or recurrent
networks is arguably much more biologically plausible than backprop. Compared
to a network of Bernoulli-logistic units trained by REINFORCE, which is argued
to be close to biological learning [1, Chapter 15], coordinated exploration
only requires adding recurrent connections between units on the same layer.

However, we only consider a one-hidden-layer ANN in this research report.
Since the coordinated exploration of both Algorithm 1 and 2 is intra-layer
instead of inter-layer, we expect that additional work is required to match
the speed of STE backprop when applied in a multi-layer ANN. It may be
possible to design a reward signal specific to each layer, with each layer
being a Boltzmann machine or a recurrent layer that receives inputs from the
previous layer. In this way, we have a vector reward signal with each reward
signal targeting a subset of neurons, and the subset of neurons communicate
with one another to generate a collective action such that a single reward is
sufficient for assigning credit to the subset; this may be closer to the
vector reward prediction errors observed in biological systems [9, 10].

Finally, further work can be done to understand how the long-term depression
(LTD) in R-STDP can benefit learning, which is not included in the learning
rules of both Algorithm 1 and 2. R-STDP states that if the pre-synaptic spike
arrives after the post-synaptic spike, the synaptic strength is decreased
instead of increased. This rule can be incorporated directly in the learning
rules of Algorithm 1 or 2, but it remains to be seen whether it can benefit
learning and what the theoretical basis is. Most studies on how R-STDP induces
learning behavior [11, 12] focus on long-term potentiation (LTP) instead of
LTD, and so how LTD can help learning is an important area for further work.

##  3 Acknowledgment

We would like to thank Andrew G. Barto, who inspired this research and
provided valuable insights and comments.

##  Appendix A Additional Experiments

###  A.1 Coordinated Exploration with Boltzmann Machines

To see how the performance of Algorithm 1 scales with $T$, the number of
sampling step, we repeat the experiments shown in Figure 3 but with $c=0.25$
and $T\in\\{1,2,4,10,15,20,25,30\\}$. The results are shown in Figure 7. We
observe that in terms of learning speed, $T\geq 15$ works well; in terms of
asymptotic performance, $T=25$ works the best. It is interesting to observe
that the asymptotic performance of $T=50$ is slightly worse than $T=25$. The
reason for this remains to be understood.

Figure 7: Learning curve of Algorithm 1 with varying $T$ and $c=0.25$. Results
are averaged over 5 independent runs, and the error bar represents the
standard deviation over the runs.

Finally, we consider how the optimal strength of recurrent connection $c$
changes with number of hidden units $N$. The average reward throughout all
episodes for Algorithm 1 with $N\in\\{8,16,32,64,96,128\\}$ and
$c\in\\{0.00,0.05,0.10,0.25,0.40,0.50,0.75,1.00\\}$ is shown in Figure 8. We
used $T=25$ here. We observe that the optimal $c$ scales inversely with $N$.
The optimal $c$ when $N\leq 16$ is around $0.50$ while the optimal $c$ when
$N>16$ is around $0.25$.

Figure 8: Average reward throughout all episodes of Algorithm 1 with varying
$N$ and $c$. Results are averaged over 5 independent runs, and the error bar
represents the standard deviation over the runs.

###  A.2 Coordinated Exploration with Recurrent Networks

To see how the performance of Algorithm 2 scales with $T$, the number of
sampling step, we repeat the experiments shown in Figure 6 but with $c=0.25$
and $T\in\\{1,2,4,8,10,20\\}$. The results are shown in Figure 9. We observe
that there are almost no differences in the learning curves for different $T$.

Figure 9: Learning curve of Algorithm 2 with varying $T$ and $c=0.25$. Results
are averaged over 5 independent runs, and the error bar represents the
standard deviation over the runs.

To see how the performance of Algorithm 2 scales with $\lambda$, the trace
decay rate, we repeat the experiments shown in Figure 6 but with $c=0.25$ and
$\lambda\in\\{0.00,0.10,0.25,0.40,0.50,0.75,0.90,0.99\\}$. The results are
shown in Figure 10. We observe that there are almost no differences in the
learning curves for $\lambda\leq 0.50$, indicating that we can apply REINFORCE
on the last time step and the performance is still the same. However,
performance deteriorates as $\lambda$ grows larger than $0.50$ due to the
increased noise.

Figure 10: Learning curve of Algorithm 2 with varying $\lambda$ and $c=0.25$.
Results are averaged over 5 independent runs, and the error bar represents the
standard deviation over the runs.

Finally, we consider how the optimal strength of recurrent connection $c$
changes with the number of hidden units $N$. The average reward throughout all
episodes for Algorithm 2 with $N\in\\{8,16,32,64,96,128\\}$ and
$c\in\\{0.00,0.05,0.10,0.25,0.40,0.50,0.75,1.00\\}$ is shown in Figure 11. We
used $T=2$ and $\lambda=0.25$ here. We observe that the optimal $c$ scales
inversely with $N$. The optimal $c$ for $N=8,16,32,64,96,128$ are respectively
$1.00,1.00,0.75,0.25,0.25,0.10$.

Figure 11: Average reward throughout all episodes of Algorithm 2 with varying
$N$ and $c$. Results are averaged over 5 independent runs, and the error bar
represents the standard deviation over the runs.

###  A.3 Miscellaneous

For completeness, we repeat Figure 4 and 5 here, but we show the average
reward throughout all 4e6 episodes instead of only the first or the last 1e6
episodes. The results are shown in Figure 12.

Figure 12: Average reward throughout all episodes for a network with varying
$N$. Results are averaged over 5 independent runs, and the error bar
represents the standard deviation over the runs.

## References

  * [1] Richard S Sutton and Andrew G Barto  “Reinforcement learning: An introduction”  MIT press, 2018
  * [2] Stephen Chung  “MAP Propagation Algorithm: Faster Learning with a Team of Reinforcement Learning Agents”  In _Advances in Neural Information Processing Systems_ 34, 2021
  * [3] Stephen Chung  “Learning by competition of self-interested reinforcement learning agents”  In _Proceedings of the AAAI Conference on Artificial Intelligence_ 36.6, 2022, pp. 6384–6393
  * [4] David H Ackley, Geoffrey E Hinton and Terrence J Sejnowski  “A learning algorithm for Boltzmann machines”  In _Cognitive science_ 9.1 Elsevier, 1985, pp. 147–169
  * [5] Wulfram Gerstner, Werner M Kistler, Richard Naud and Liam Paninski  “Neuronal dynamics: From single neurons to networks and models of cognition”  Cambridge University Press, 2014
  * [6] Yoshua Bengio, Nicholas Léonard and Aaron Courville  “Estimating or propagating gradients through stochastic neurons for conditional computation”  In _arXiv preprint arXiv:1308.3432_ , 2013
  * [7] Ronald J Williams  “Simple statistical gradient-following algorithms for connectionist reinforcement learning”  In _Machine learning_ 8.3-4 Springer, 1992, pp. 229–256
  * [8] Diederik P Kingma and Jimmy Ba  “Adam: A method for stochastic optimization”  In _arXiv preprint arXiv:1412.6980_ , 2014
  * [9] Nathan F Parker et al.  “Reward and choice encoding in terminals of midbrain dopamine neurons depends on striatal target”  In _Nature neuroscience_ 19.6 Nature Publishing Group, 2016, pp. 845–854
  * [10] Rachel S Lee, Ben Engelhard, Ilana B Witten and Nathaniel D Daw  “A vector reward prediction error model explains dopaminergic heterogeneity”  In _bioRxiv_ Cold Spring Harbor Laboratory, 2022
  * [11] Răzvan V Florian  “Reinforcement learning through modulation of spike-timing-dependent synaptic plasticity”  In _Neural Computation_ 19.6 MIT Press, 2007, pp. 1468–1502
  * [12] Nicolas Frémaux, Henning Sprekeler and Wulfram Gerstner  “Reinforcement learning using a continuous time actor-critic framework with spiking neurons”  In _PLoS computational biology_ 9.4 Public Library of Science, 2013


# Local Field Potential Journey into the Basal Ganglia

Eitan E. Asher1, Maya Slovik2,3, Rae Mitelman2,3, Hagai Bergman2,3, Shlomo
Havlin1, Shay Moshel1,2,3,4 1Department of Physics, Bar-Ilan University, Ramat
Gan, Israel
2Department of Medical Neurobiology, The Hebrew University-Hadassah Medical
School, Jerusalem, Israel
3 The Edmond and Lily Safra Center for Brain Sciences, The Hebrew University,
Jerusalem, Israel 4 Nuclear Research Center Negev, Beer-Sheva, Israel

###### Abstract

Local Field potential (LFP) in the basal ganglia (BG) nuclei in the brain has
attracted much research and clinical interest. However, the origin of this
signal is still under debate throughout the last decades. The question is
whether it is a local sub threshold phenomena, a synaptic inputs to neurons or
it is a flow of electrical signals merged as volume conduction which are
generated from simultaneous firing neurons in the cerebral cortex and obeys
the Maxwell equations. In this study, we recorded in monkey brain
simultaneously LFP’s from the cerebral cortex, in the frontal lobe and primary
motor cortex (M1) and in sites in all BG nuclei: the striatum, globus pallidus
and subthalamic nucleus. All the records were taken from human primate model
(vervet monkey), during spontaneous activity. Developing and applying a novel
method to identify significant cross correlations (potential links) while
removing ”spurious” correlations, we found a tool that may discriminate
between the two major phenomena of synaptic inputs (as we define as
information flow) and volume conduction. We find mainly two major paths flows
of field potential, that propagate with two different time delays, from
primary motor cortex and from the frontal cortex. Our results indicate that
the two paths flows may represent the two mechanisms of volume conduction and
information flow.

Local Field Potential, Links

##  I Introduction

Broad band field potential becomes an important and powerful tool to
understand many neuronal functions, abnormalities and clinical phenomena.
Behavioral studies use local field potential (LFP) activity response to
describe cognitive occurrences Buzsáki et al. (2012); Wimmer et al. (2016);
Widge et al. (2019) and sensory response Henrie and Shapley (2005); Bosman et
al. (2012). Brain machine interface tools have been developed based on
computation of local field potential Jackson and Hall (2016) and
electroencephalogram (EEG) signals Müller et al. (2008). Subthalamic nucleus
LFPs has been found to be coupled with EEGs related to cortical somatosensory
areas Marsden et al. (2001),and neural discharge rates found to be coupled
with LFPs signals in the cerebral cortex during voluntary movements Donoghue
et al. (1998). Brain pathology diseases such as Parkinson disease (PD) and
schizophrenia, cause to the increase of brain abnormal oscillations in
neuronal signals, as single-unit, multi-unit activity and LFPs in the brain,
such as, in the Basal Ganglia (BG) and cerebral cortex Rodriguez-Oroz et al.
(2001); Moshel et al. (2013); Slovik et al. (2017). LFPs span the frequency
ranges of 1-80 Hz or 1-400 Hz, if one includes the high gamma peaks reported
at 65-90 Hz and 250-350 Hz Danish et al. (2007). However, Yuval-Greenberg et
al. (2008) shows possible confounding factors in the high frequency regime of
LFP, whereas spikes have their maximal power around 1000 Hz. Thus, although
LFP oscillations have been thought to imply spike synchronization Brown and
Williams (2005); Hammond et al. (2007); Özkurt et al. (2011); Moshel et al.
(2013), studies claim that they are more likely to represent sub-threshold
phenomena as synaptic inputs Belitski et al. (2010); Buzsáki et al. (2012)
which is probably correlated with spike activity. Broadband field potentials
have been also found to be correlated with spontaneous cerebral low-frequency
blood oxygen level-dependent (BOLD) fluctuations Pan et al. (2011); Hermes et
al. (2017). Hermes et al described that, the two measures, the BOLD amplitude
and electrocortiogram (ECoG) broadband power are correlated in cerebral cortex
(V1, V2 and V3), where the BOLD amplitude and alpha power (8–13 Hz) are
negatively correlated, and the BOLD amplitude and narrow–band gamma power
(30–80 Hz) are uncorrelated. They pointed out Hermes et al. (2017) that the
two measures, nevertheless, provide complementary information about human
brain activity, and they infer that features of the field potential that are
uncorrelated with BOLD arise largely from changes in neuronal synchrony,
rather than, level of direct neuronal activity. Although it is found that
broadband field potential plays a key role for analyzing brain behavior and it
is easy to record, the origin of the local field potential (LFP) remains under
debate. Many studies suggest that this signal carries a valuable neuronal
information as the local synaptic inputs. In spite of these, there are studies
that describe these signals as propagative electromagnetic waves behave as a
volume conductance, obeying the Maxwell equations. An intermediate studies
suggest that these signals have mixed characters of both volume conduction and
a summation of the extracellular fields Buzsáki et al. (2012). Parabucki et al
Parabucki and Lampl (2017) described a strong LFP signals coupling between the
cerebral cortex, during whisker stimulation in a mice, with a dissected part
of the olfactory bulb. They recorded simultaneously the cerebral cortex and
the dissected part, and describe that the LFP “passed” from the cortex through
the dissected part as a volume conducted signal. However, in spite of these
remarkable achievements, to date we still do not have a clear picture for
describing the activity of broadband field potentials, whether it is a volume
conducted electromagnetic waves that spreads inside the deep layers in the
brain or they carry a valuable information as synaptic inputs. In order to
explore this issue, we map LFP signals from the cerebral cortex trough the
basal ganglia nuclei, where we recorded simultaneously the signals from the
cerebral cortex and signals in the deep layers of the basal ganglia (BG)
nuclei, of human primate model during wake. Using micro-electrodes, we
recorded the full electrophysiological signal (row-data) in the range of [1,
6000] Hz, with 1-7 electrodes which were fixed in the 5th cortex layer
(pyramidal cells) for all recording session, see Fig.1(a). One electrode that
has been moved towards to the BG nuclei and has been simultaneously recorded
every 200 $\mu$m, for a duration of 120 seconds, since the entrance to the
cerebral cortex (frontal lobe) and passed through: cortex, striatum, globus
pallidus (external - GPe, internal - GPi), internal capsule and the
subthalamic nucleus. To get LFP signal we filtered the full signal (raw data)
to 0.1-80Hz and analyzed the relationship between neuronal activity in the
cerebral cortex and the deep BG brain layers.

##  II Methods

###  II.1 Animals, surgery and MRI

We conducted the experiments on one African Green monkey (Cercopithecus
aethiops, females, weighing between 3.3 and 4.4 Kg). All procedures were
conducted in accordance with the Hebrew University guidelines for animal care
and the National Institute of Health Guide for the Care and Use of Laboratory
Animals. The Hebrew University is an AAALAC (Association for Assessment and
Accreditation of Laboratory Animal Care) approved institute. The monkey was
trained to sit quietly in a primate chair. After training, the animal
underwent a surgical implantation of a Cylux MRI compatible 27*27 mm recording
chamber (Alpha-Omega, Israel) and head holder (Crist Instruments, MD). The
location of the chamber was determined by a primate stereotaxic atlas
(Contreras et al. 1981; Martin et al. 2000) to be above the right primary
motor cortex (M1), glubos pallidus (Gp) and the subthalamic nucleus (Stn). A
few days after surgery, an MRI scan was performed to determine its exact
stereotaxic location.

###  II.2 Recording and data acquisition

During recording sessions, the animals were seated in a primate chair with
their head and hands restrained, but were free to move their trunk and legs.
All recording were done during alert monkey’ spontaneous behavior. The
verification of recording location was carried out separately for the target
structures, the cerebral cortex and the deep nuclei of the basal ganglia (BG)
complex: striatum, globus pallidus and the subthalamic nucleus. Extra cellular
data domain was amplified with a gain of 5K, band-pass filtered with a 1-6000
Hz four-pole hardware Butterworth filter and continuously sampled at 25 kHz.

Figure 1: Hemi-coronal illustration from Macaque AC atlas - 6 mm and
electrodes setup, raw data and power spectrum densities. (a) We place 7
electrodes in the cerebral cortex, fixed on the 5th layer of pyramidal cells.
Electrodes 1-4 (E1-4) were permanently located on the right–hemisphere,
primary M1-motor-cortex, and electrodes 5-7 (E5-7) on the right hemisphere -
frontal cortex. One electrode aimed to the basal ganglia, moving to the
deepest layers of the basal ganglia passing through the recording sites in:
the cerebral cortex (Ctx), striatum (Str), globus pallidus (Gp) internal
capsule (Ic) and subthalamic nucleus (Stn). We start the recording in the
entrance of the nucleus surface as been detected according to a typical
electrophysiological activity. Each site has been recorded for time duration
of 120 seconds. After completing the recording, we move in a step of 200
micrometer deeper, as illustrated by the arrows. (labels in cartoon: A to F
mark the recorded nucleus). (b) Examples of 5 seconds records of raw local
field potential data as recorded from the nucleus that were analyzed, from
shallow to deep: Ctx, Str, Gp, Ic and Stn, and (c) their power spectral
density, after applying a band-pass filter to remove the 50Hz band.

The data has been stored by a data acquisition system (Alpha-Map, Alpha-Omega,
Israel) for offline analysis. The animals’ behavior has been monitored by a
computerized video surveillance system (monkeys F and Su: GV-650, GeoVision,
Taiwan; monkeys F and Su: GV-2400, Geo Vision, Taiwan) connected to four
infrared digital cameras. We used two arrays of four glass-coated tungsten
microelectrodes (impedance measured at 1 kHz, range 0.2-0.6 MOhm), confined
within a cylindrical metal guide (1.36 mm internal diameter). The electrodes
were positioned (E1-4 electrodes) in the area of the primary motor cortex and
the other E5-7 electrodes were fixed at frontal cortex 5 th layer and aimed to
the basal ganglia nuclei , see Fig.1(a). This was done using the double
Microdriving Terminal and Electrode Positioning System (Double MT and EPS,
Alpha-Omega, Israel), that allows independent movement of each electrode.
During the recording procedure, we fixed E1-7 electrodes at the 5th layer of
the cerebral cortex (M1 and the frontal cortex above the BG) and a moving 8th
electrode, has been advanced towards the BG nuclei (Figure 1A). We started
simultaneous recordings of the Ctx electrodes (E1-7) when the moving electrode
entered the cerebral cortex. We recorded at each location for time duration of
two minutes, after waiting several seconds for stabilization. When we
completed the specific site recording; we moved further the traveling
electrode by 200 $\mu$m to the next deeper location. The verification of
recording location was carried out separately for the target structures, using
their typical electro-physiological signatures (Moshel et al., 2013). In our
data analysis we utilized Python 3.6 and custom-made MATLAB19a (R2019a)
routines. Our database contains 13 trajectories from the basal ganglia nuclei.
Data records are summed in Table 1.

Location | E1-7-Ctx | E1-7-Str | E1-7-Gp | E1-7-Ic | E1-7-Stn
---|---|---|---|---|---
Sites | 210 | 219 | 91 | 20 | 18
Pairs | 1318 | 1331 | 561 | 140 | 126
Table 1: Data of the study. Number of recorded sites of each nucleus, and
number of total cortex electrode-deep electrode pairs. The deep electrode was
moved in 200$\mu$m steps from the cortex towards the BG. Each location that
has been visited is a site. At each site, the signal that was recorded from
the deep electrode has been analyzed using cross correlation against its
matching time record of each of the E1-7 electrodes, resulting in the number
of total pairs, given in the Table.

###  II.3 Cross-correlation and surrogate analysis to identify significant
interactions

In order to quantify the interaction represented by similarity between two
signals in both sites X and Y. (One site is always the site from the deep
electrode, and the other one is taken from the electrodes that are located on
the cortex (electrodes E1-7), see Fig.1. To this end, we analyze similarity,
using cross correlation, between the of two simultaneous signals ${x}$ and
${y}$ each of length N. We divide both signals ${x}$ and ${y}$ into
$N_{L}$-overlapping segments $\nu$ of equal length L=2.5 s. We choose an
overlap of L/4=0.625 s. Thus, for each 120 second record we get $N_{L}=63$
time windows. Before the analysis, the signal in each segment $\nu$ is
subtracted by its mean, in order to remove constant trends in the data. This
normalization assures that the coupling between the signals ${x}$ and ${y}$ is
not affected by their relative amplitudes, but mainly on variations from
average.. Next, we calculate the cross correlation function,
$R^{\nu}_{xy}(\tau)=\sum_{\tau=-n/2}^{n/2}x^{\nu}(n+\tau)*y^{\nu}(n)$ within
each segment $\nu=1,...,N_{L}$. We then normalize the sequence so that the
auto-correlations at zero lag equals 1:
$\hat{R}(\tau)=\frac{1}{\sqrt{\hat{R}_{xx}(0)\hat{R}_{yy}(0)}}\hat{R}_{xy}(\tau)$
for each window $\nu$. In order to distinguish between significant and non-
significant interactions between signals ${x}$ and ${y}$, we study how their
cross-correlation $R$ decays when shifting the signals against each other. The
more it decays the more significant are the links. Panels (e) and (f) in Fig.
(2) show examples of R versus time shift $\tau$ and suggest that more
synchronized signals (with higher R values) have a marked decay of R($\tau$)
(panel (e)) that may not be seen for less synchronized signals (panel (f)). To
better quantify this observation, we define a significance value $w$ (standard
score) that normalizes the maximum cross correlation index $R_{max}$ (detected
at a particular time shift $\tau^{*}$) by the “background noise” characterized
by the mean and standard deviation of R($\tau$). Colloquially, $w$ gives an
estimate of how much $R_{max}$ “stands out” from the background. Specifically,
how many standard deviations the value of the peak is farther from the noise.
That is,

$w^{i,j}=\frac{(|R_{max}|)-\langle R\rangle)}{\sigma(R)}$ (1)

where $\langle R(\tau)\rangle$ ($\sigma(R(\tau))$) is the mean (standard
deviation) of $R(\tau)$. From Eq. (1) one can see that if $R_{max}$ is much
larger than the background, $w$ is higher, indicating a more significant
coupling between signals $x$ and $y$. On the other hand, if $R_{max}$ does not
stand out from the $R(\tau)$ background, $w$ is lower and implies non
significant $x–y$ interaction. Figure (2)g) shows results obtained from 10,000
pairs of signals, where each pair of signals is either taken at the same time
(simultaneously, real” data) or at different times (“surrogate” data). It can
be seen from the figure that high $w$ values for the real data are usually
detected for small shifts $\tau^{*}$ (of the order of few tens of
milliseconds), which is consistent with the fact that brain waves generally
propagate rather quickly Zhang et al. (2018); Sato et al. (2012). Combining
the two conditions, we consider $x–y$ interactions as significant only if
$\tau^{\ast}\in[-0.05,0.05]$ seconds and $w>4.5$.

Figure 2: Correlations and surrogate analysis for identifying significant
interactions. (a-f) Demonstration of two pairs of signals (blue curves in
(a)-(c) and (b)-(d)) obtained simultaneously from different nuclei after
applying a [0-80Hz] low-pass-filter to the raw LFP data. (e) and (f) The cross
correlation index $R$ as a function of the shift $\tau$ between the signals
(a) vs. (c) and (b) vs. (d), respectively. The signals from (a) and (c) yield
a marked maximum $R$ at shift $\tau^{\ast}=\left.\tau\right|_{R(\tau)\equiv
R_{\rm{max}}}=0$, and $R(\tau)$ decays rapidly for $|\tau|>0$. (b) and (d)
demonstrate uncorrelated signals, since, $R(\tau)$ in (f) shows fluctuating
behavior without clear decay, i.e, uncorrelated signals. A significance value
$w$ characterizes $R(\tau)$ by normalizing $R_{\rm{max}}$ by the mean and
standard deviation of $R(\tau)$ (Eq. 1). Correspondingly, we obtain a higher
$w$ value for panel (e) ($w=5.2$) than for panel (f) ($w=.67$). (g)
distribution of $w$ for all data, taken from electrodes and E5-7, compared to
surrogate signals. For the real data, we show $w$ for all time windows that
were analyzed. Surrogate signals are signals that were randomly chosen from
different electrodes, different nucleus, and also at different times. Note
that the E5-7 (frontal-lobe) electrodes, the distribution of real $w$ ’s is
higher, Compared to those of the ones were fixed at cortex M1 , electrodes
E1-4, which means that Lfps from frontal lobe are more correlated. Surrogate
signals yield $w$ values distribution that is significantly lower than $w$’s
of real signals; maximum $w$ value for surrogate signals is around 4. Panel
(h) shows that highest $w$ values are observed for $\tau^{\ast}\approx 0$. In
this scatter plot we show 5000 samples of $w$ vs. $\tau^{\ast}$ for real data
(blue circles) and surrogate data (red dots). Real signals are taken at the
same time, comparing deep and Ctx signals, whereas surrogate pairs were chosen
randomly from different times. Clearly, higher $w$ values are obtained for
real signals for $\tau^{\ast}\approx 0$. In contrast, surrogate analysis does
not lead to higher $w$ values around $\tau^{\ast}\approx 0$ and shows a
uniform $w$ vs. $\tau^{\ast}$ distribution.

##  III Results

Correlations, weights and $\tau^{*}$ distributions. We find that there exist
strong correlations between Ctx-Ctx, Ctx-Str, Ctx-Gp, Ctx-Ic and Ctx-Stn, see
Fig. 11. We define the time lag $\tau^{\ast}$ at which $C(\tau)$ is maximal,
as the time delay of pairs $x$ and $y$. Whether $\tau$ is positive or
negative, determines the location in which the signal was detected first. We
find that when looking on real links in which $w>4.5$, the
$\langle\tau^{\ast}\rangle$ is around factor two larger for the electrodes
located at primary motor cortex - M1 (E1-4) than the ones located at frontal
cortex E5-7, see Fig. 4 (a). Additionally, when plotting the distributions of
$\tau^{\ast}$ for all possible interactions, the distribution of E1-4
electrodes is much wider than E5-7.

Link lifetime. Another question of interest is what is the time duration or
”lifetime” of a link? As before, we divide each recording time period at each
site (120 s) into 2.5 s overlapping windows, with a fraction of 0.25 window
length overlap. Then, by calculating the cross-correlation of the signal of a
given site with its simultaneous time signal in the cortex electrode, we find
the link strength $w$ by applying Eq.1. If $w>4.5$, we regard the two signals
as linked. Then, we find the time duration (how many successive windows are
correlated) of those links, and plot their distribution for all the nuclei
interactions that were probed. The results are presented in Fig.5. We find
that, as the interaction involves towards deeper nucleus, there is a clear
decrease in the link length.

Figure 3:  Cross-correlation examples. a. A cross-correlation result of two
signals, one is taken from the deep electrode and the second is from electrode
number 7, which is located on the frontal lobe. The $\tau^{\ast}$ value
$<0.02$ which corresponds to a first mode interaction. In b, same as a, where
the second signal is obtained from electrode number 2 which records the M1
motor cortex. Here, $\tau^{\ast}=0.03$ which corresponds to a second mode
interaction.

Modes analysis. Next, we observed a bi-modal distribution of $\tau^{\ast}$
(See Fig. 4, b), supporting the hypothesis of the existence of two transfer
mechanisms. The first one, dominated by $\tau^{\ast}$ centered around zero
time lag, indicating a volume conduction, and a second mode, with larger
$\tau^{\ast}$ (typically 30 milliseconds) could be related to information
transfer through synaptic inputs. As shown above, every significant link is
composed of several consecutive sub-links of 2.5 seconds each. To test the
above hypothesis of the existence of two modes we will test next if all or at
least most of sub-links belong to the same mode or if we have a random mixture
of modes. If the links are mainly unique that is almost all sub-links in a
specific link have the same $\tau^{\ast}$, as indeed found here, it supports
the hypothesis of two modes. To this end, we examined the links sequences
compositions, i.e., the sub-links, by dividing the $\tau^{\ast}$ range in Fig.
4 b into two regimes by the interception points of two fitted Gaussians, which
we denote by ‘1’ if $\tau^{\ast}$ of the window belongs to the first regime
$\tau^{\ast}$, and ‘2’ if it falls in the second regime $\tau^{\ast}$.

Then, for each sequence, we find the sequence composition of sub-links. For
example, a sequence of 3 sub-links, may result in the composition ‘1,1,2’,
meaning that the first and second windows had $\tau^{\ast}$’s from the first
mode, and third window had $\tau^{\ast}$ in the second mode regime, as
illustrated on Fig. 6. Interestingly, we find that for most of the sequences,
the whole sequence was composed by either consecutive ‘1’ or consecutive ‘2’.
Meaning, that every interaction is of one type, and dominated by one transfer
mechanism. We assume that the first mode, which consists of $\tau^{\ast}\sim
0$, are interactions that occur due to volume conduction process, and that the
second mode of $\tau^{\ast}\sim 0.02$ represent information flow. It can be
seen that only a very small fraction of the sequences are composed of mixed
modes of 1’s and 2’s, (e.g., ‘1,2,1’). The results are summarized in Fig. 6.

Figure 4: Basal ganglia- cerebral cortex, average time delay results.  (a)
Average $|\tau^{\ast}|$’s of nuclei, of interactions with $\tau>0$. We show
the mean time delays of the two signals of the electrodes sets, i.e., between
E1-4, E5-7, and the moving electrode. $\tau>0$ corresponds to the case in
which the signal initiated (or earlier) at the cortex. This consist $80\%$ of
the real links. Note that the electrodes located at M1 (E1-4) have systematic,
for all nuclei larger time delays, than the frontal lobe electrodes (E5-7)
with a similar profile. (b) The distribution of $\tau^{\ast}$ for all cerebral
cortex-BG nuclei Note that the distributions are bi-modal, suggesting two
types of links. In light green, we present results of surrogate Lancaster et
al. (2018) data, where the deep electrode is taken form the nucleus, and the
cortical electrodes are taken from different time of recording. Surrogate
$\tau^{\ast}$ shows a uniform distribution without particular clustering
behavior at specific $\tau^{\ast}$ as expected (light green). These results
suggests that the mechanisms of volume conductance and information transfer
occur. The number of time windows that were analyzed per nuclei appear on the
top left box of each nuclei.  Figure 5:  Distribution of links duration of
different nuclei.  Left and right columns show PDFs of links duration of E1-4
and E5-7 electrode sets respectively. Each site record (120 sec.), is divided
in to 2.5 seconds long overlapping windows, with $25\%$ overlap (0.625
seconds). Then, we analyze the cross-correlation within each time window
segment between signals in the deep electrode and its matching time record of
the electrode in the cortex, and $w$ is calculated. If $w>4.5$, we assume that
the link is real. We then find the duration of the consecutive real links that
were formed with the cortex electrodes at all sites of record. The results
presented in log-log scale, suggest a power law distribution behavior
(exponent $\gamma$ is shown), and declination of real links duration as the
nuclei is deeper, is visible. Note also, that for most nuclei links lengths
are longer for the E5-7 frontal lobe electrode set compared to E1-4 motor M1
set. This suggests that distance affects LFP signals. Exponents show
consistent increase with depth, except Ic, which consists of white matter in
which signals might propagate more rapidly.  Figure 6:  Illustration of the
first, the second and the mixed modes sequences composition. This figure
comprise of three rows. The bottom row represents a 1st mode sequence of
length $n$. Note that for this sequence $\tau^{\ast}=0$ for all windows. A
sequence is $n$ time-windows of which each time window had $w>4$. Each
sequence of length $n$, consists of $n$ time windows that were cross-
correlated. We calculate the $\tau^{\ast}$ of each of the $n$ time windows,
and ask what is the composition of those $\tau^{\ast}$’s, i.e, to what mode
this window belongs.  Figure 7: Sequences modes composition, $w>4.5$.  (a)
analysis of compositions of links obtained between frontal cortex electrodes
and the deep electrode. (b): same as top, for the M1 electrodes. For each link
of time length L, we analyze its sub-links $\tau\ast$ values, and test whether
they are in the ranges of the first or second modes. We find, that for most of
the sequences, the sequence belongs to the the 1st or 2nd modes exclusively,
while few of the sequences contains a mix of the both modes. If the whole
sequence contains only $\tau\ast$ values from the first (second) mode, we say
that this sequence is a ”first (second) mode” sequence. If the sequence
contains combinations of first and second modes, we say it is a ”mixed”
sequence. Note that the frontal cortex links are dominated by first mode
sequences $\langle\tau\rangle\backsimeq 0$, while second mode sequences
$\langle\tau\rangle\backsimeq 0.02$ appear much more in the M1 links. Also
note, that for the internal capsule (Ic) which consists of white matter, most
links are from first mode in both frontal cortex and M1.

##  IV Summary and discussion

In this work, We tested the open question whether the LFP signals is a volume
conductance phenomena or a synaptic input which represents information
transfer. To adress this question we develop a significant based cross
correlation method and define lifetime of links between different locations of
the brain. To adress this question we develop a significant based cross
correlation method and define a continuous time delayed link patterns. We
performed simultaneous recordings of local field potentials (LFP) in the range
of [0.1, 6000] Hz by 7 micro-electrodes which have been fixed on the 5th layer
in the cerebral cortex and we used one micro-electrode to record sites while
moving towards the BG nuclei. Each site was $200\mu m$ deeper than the
previous one, and was aimed towards the BG. This procedure was done since the
entering of the deep electrode to the frontal lobe cortex until it reached to
the subthalamic nucleus. See illustration in Fig. 1(a). After data were
collected and digitized, power spectrum densities were computed using Welch’s
power spectral density estimate method (Fig. 1(c)). Since the frequency range
that was recorded includes the local field potential and the single and multi-
unit activity, we apply low-pass filter to the signals in the range of [0.1,
80] Hz in order to get the field potentials. We split the 120 second duration
of the recording at each site to 2.5 second windows with $0.5$ seconds overlap
between each window, and compute the cross correlation between each recording
signal of the deep electrode with its matching time window from the seven
cerebral cortex. Next, we developed a significance weight method to
distinguish between spurious and real links. By applying this to a surrogate
control data, we found a threshold score that enables us to assign for each
link a weight which corresponds to its strength, and we compared the weights
distribution to that of the surrogate signals in order to identify the real
links. After removing the spurious interactions, we also find the time delay
$\tau^{\ast}$ for each real link using the peaks in the shifting cross
correlations. When plotting the average $\tau^{\ast}$ for both M1 and cortex
frontal - lobe for each nuclei, we see (Fig. 4) that frontal-cortex areas
yield shorter $\tau^{\ast}$. Moreover, frontal-cortex area also produce more
significant interactions, which can be seen in their weight scores (Fig.
2(g)).
When we compared real links duration from M1 and the cortex, we observed that
as the recorded nuclei is deeper, links durations become shorter, and links of
the deep electrode-cortex are longer than their corresponding deep
electrode-M1, see Fig. 5. Furthermore, when looking at the $\tau^{\ast}$
distributions, we surprisingly observe two main modes, see Fig. 4(b). One,
characterized by a short time delays ($\tau^{\ast}\sim 0$), and the second,
with longer time delays ($\tau^{\ast}\sim 0.02$). We hypothesize, that these
two modes are related to two different mechanisms of interaction, volume
conductance (short $\tau^{\ast}$, 1st mode) and threshold phenomena as
synaptic inputs - an information transfer (longer $\tau^{\ast}$, 2nd mode).

Our next question was to test whether a long term interaction, e.g. link, is
mainly characterized by a $\tau^{\ast}$ sequence that belongs to one of the
two modes. Interestingly, we find that most of the links consist of mainly 1st
or 2nd mode only, and very few are a mix of the two, see Fig. 7. Frontal
cortex - deep electrode links comprise mostly of 1st mode sequences, while 2nd
mode sequences are more prevalent in the M1 - deep electrode links.

The significance of our method and results is that it may provide a framework
to distinguish between volume conduction interactions, and real information
flow, which currently an unsolved question. Future work may use tools from
information theory and test whether links that belong to the second mode,
conveys more information compared to those of volume conductance.

##  V Acknowledgments

We thank..

## References

  * Buzsáki et al. [2012] György Buzsáki, Costas A Anastassiou, and Christof Koch.  The origin of extracellular fields and currents—eeg, ecog, lfp and spikes.  _Nature reviews neuroscience_ , 13(6):407–420, 2012.
  * Wimmer et al. [2016] Klaus Wimmer, Marc Ramon, Tatiana Pasternak, and Albert Compte.  Transitions between multiband oscillatory patterns characterize memory-guided perceptual decisions in prefrontal circuits.  _Journal of Neuroscience_ , 36(2):489–505, 2016\.
  * Widge et al. [2019] Alik S Widge, Sarah R Heilbronner, and Benjamin Y Hayden.  Prefrontal cortex and cognitive control: new insights from human electrophysiology.  _F1000Research_ , 8, 2019.
  * Henrie and Shapley [2005] J Andrew Henrie and Robert Shapley.  Lfp power spectra in v1 cortex: the graded effect of stimulus contrast.  _Journal of neurophysiology_ , 94(1):479–490, 2005.
  * Bosman et al. [2012] Conrado A Bosman, Jan-Mathijs Schoffelen, Nicolas Brunet, Robert Oostenveld, Andre M Bastos, Thilo Womelsdorf, Birthe Rubehn, Thomas Stieglitz, Peter De Weerd, and Pascal Fries.  Attentional stimulus selection through selective synchronization between monkey visual areas.  _Neuron_ , 75(5):875–888, 2012.
  * Jackson and Hall [2016] Andrew Jackson and Thomas M Hall.  Decoding local field potentials for neural interfaces.  _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , 25(10):1705–1714, 2016.
  * Müller et al. [2008] Klaus-Robert Müller, Michael Tangermann, Guido Dornhege, Matthias Krauledat, Gabriel Curio, and Benjamin Blankertz.  Machine learning for real-time single-trial eeg-analysis: from brain–computer interfacing to mental state monitoring.  _Journal of neuroscience methods_ , 167(1):82–90, 2008.
  * Marsden et al. [2001] JF Marsden, P Limousin-Dowsey, P Ashby, P Pollak, and P Brown.  Subthalamic nucleus, sensorimotor cortex and muscle interrelationships in parkinson’s disease.  _Brain_ , 124(2):378–388, 2001.
  * Donoghue et al. [1998] John P Donoghue, Jerome N Sanes, Nicholas G Hatsopoulos, and Gyongi Gaal.  Neural discharge and local field potential oscillations in primate motor cortex during voluntary movements.  _Journal of neurophysiology_ , 79(1):159–173, 1998.
  * Rodriguez-Oroz et al. [2001] Maria C Rodriguez-Oroz, Manuel Rodriguez, Jorge Guridi, Klaus Mewes, Vijay Chockkman, Jerrold Vitek, Mahlon R DeLong, and Jose A Obeso.  The subthalamic nucleus in parkinson’s disease: somatotopic organization and physiological characteristics.  _Brain_ , 124(9):1777–1790, 2001.
  * Moshel et al. [2013] Shay Moshel, Reuben Ruby Shamir, Aeyal Raz, Fernando Ramirez de Noriega, Renana Eitan, Hagai Bergman, and Zvi Israel.  Subthalamic nucleus long-range synchronization—an independent hallmark of human parkinson’s disease.  _Frontiers in systems neuroscience_ , 7:79, 2013.
  * Slovik et al. [2017] Maya Slovik, Boris Rosin, Shay Moshel, Rea Mitelman, Eitan Schechtman, Renana Eitan, Aeyal Raz, and Hagai Bergman.  Ketamine induced converged synchronous gamma oscillations in the cortico-basal ganglia network of nonhuman primates.  _Journal of neurophysiology_ , 118(2):917–931, 2017.
  * Danish et al. [2007] SF Danish, JT Moyer, LH Finkel, GH Baltuch, JL Jaggi, A Priori, and G Foffani.  High-frequency oscillations (¿ 200 hz) in the human non-parkinsonian subthalamic nucleus.  _Brain research bulletin_ , 74(1-3):84–90, 2007\.
  * Yuval-Greenberg et al. [2008] Shlomit Yuval-Greenberg, Orr Tomer, Alon S Keren, Israel Nelken, and Leon Y Deouell.  Transient induced gamma-band response in eeg as a manifestation of miniature saccades.  _Neuron_ , 58(3):429–441, 2008.
  * Brown and Williams [2005] Peter Brown and David Williams.  Basal ganglia local field potential activity: character and functional significance in the human.  _Clinical neurophysiology_ , 116(11):2510–2519, 2005.
  * Hammond et al. [2007] Constance Hammond, Hagai Bergman, and Peter Brown.  Pathological synchronization in parkinson’s disease: networks, models and treatments.  _Trends in neurosciences_ , 30(7):357–364, 2007\.
  * Özkurt et al. [2011] Tolga Esat Özkurt, Markus Butz, Melanie Homburger, Saskia Elben, Jan Vesper, Lars Wojtecki, and Alfons Schnitzler.  High frequency oscillations in the subthalamic nucleus: a neurophysiological marker of the motor state in parkinson’s disease.  _Experimental neurology_ , 229(2):324–331, 2011\.
  * Belitski et al. [2010] Andrei Belitski, Stefano Panzeri, Cesare Magri, Nikos K Logothetis, and Christoph Kayser.  Sensory information in local field potentials and spikes from visual and auditory cortices: time scales and frequency bands.  _Journal of computational neuroscience_ , 29(3):533–545, 2010.
  * Pan et al. [2011] Wen-Ju Pan, Garth Thompson, Matthew Magnuson, Waqas Majeed, Dieter Jaeger, and Shella Keilholz.  Broadband local field potentials correlate with spontaneous fluctuations in functional magnetic resonance imaging signals in the rat somatosensory cortex under isoflurane anesthesia.  _Brain connectivity_ , 1(2):119–131, 2011.
  * Hermes et al. [2017] Dora Hermes, Mai Nguyen, and Jonathan Winawer.  Neuronal synchrony and the relation between the blood-oxygen-level dependent response and the local field potential.  _PLoS biology_ , 15(7):e2001461, 2017.
  * Parabucki and Lampl [2017] Ana Parabucki and Ilan Lampl.  Volume conduction coupling of whisker-evoked cortical lfp in the mouse olfactory bulb.  _Cell reports_ , 21(4):919–925, 2017.
  * Zhang et al. [2018] Honghui Zhang, Andrew J Watrous, Ansh Patel, and Joshua Jacobs.  Theta and alpha oscillations are traveling waves in the human neocortex.  _Neuron_ , 98(6):1269–1281, 2018.
  * Sato et al. [2012] Tatsuo K Sato, Ian Nauhaus, and Matteo Carandini.  Traveling waves in visual cortex.  _Neuron_ , 75(2):218–229, 2012.
  * Lancaster et al. [2018] Gemma Lancaster, Dmytro Iatsenko, Aleksandra Pidde, Valentina Ticcinelli, and Aneta Stefanovska.  Surrogate data for hypothesis testing of physical systems.  _Physics Reports_ , 748:1–60, 2018.

##  VI Supporting information

add average tau for $\tau<0$. add results of 5 seconds window.
add correlations of the closer electrodes (in the paper only the far are
shown)
add ¡tau¿ figure for different w’s.
add links lengths figure for different w’s.

Figure 8: Note that Ctx-Ctx interactions (blue) yield higher $w$’s. Deeper
Stn-Ctx (green) interaction yield the lowest $w$ distribution, while Str-Gp-Ic
$w$ distributions increase gradually in between.  Figure 9: Average time delay
results for nuclei.  (a) Average $|\tau^{*}|$’s of nuclei, of interactions in
which $\tau<0$. We show the time delays of the two electrodes sets, M1 and
E5-7. $\tau<0$ corresponds to the case in which the interaction initiated in
the basal ganglia. Similar to the results that are shown in Fig.4(b) M1
electrodes have a larger time delays, than the E5-7 with a similar profile.
This holds for all nuclei.  Figure 10: Sequences modes composition, $w>4$.
(a) analysis of compositions of links obtained between M1 electrodes and the
deep electrode. (b): same as top, for the frontal cortex electrodes. For each
link of time length L, we analyze its sub-links $\tau\ast$ values, and test
whether they are in the ranges of the first or second modes. We find, that for
most of the sequences, the sequence belongs to the the 1st or 2nd modes
exclusively, while few of the sequences contains a mix of the both modes. If
the whole sequence contains only $\tau\ast$ values from the first (second)
mode, we say that this sequence is a ”first (second) mode” sequence. If the
sequence contains combinations of first and second modes, we say it is a
”mixed” sequence. Note that the frontal cortex links are dominated by first
mode sequences, while second mode sequences appear much more in the M1 links.
Also note, that for the internal capsule (Ic) which consists of white matter,
most links are from first mode in both frontal cortex and M1.  Figure 11:
Cross-correlations. The cross-correlations of the E5-7 (on the cortex)
electrodes with the deep electrode, while being in the five nuclei. (a-e) show
the cross-correlations between E5-7-Ctx, E5-7-Str, E5-7-Gp, E5-7-Ic and
E5-7-Stn respectively. Here we show only cross-correlations with $w>4.5$ and
$\tau\ast<0.05$ ms, of the time windows that were analyzed. In blue is the
mean, and the shaded area represents the standard deviation. Shown here are
results of E5-7 electrodes interaction. (See SI for results of E1-4
electrodes.)