# Supplementary Information: Understanding the computation of time using
neural network models

Zedong Bi  Changsong Zhou

##  S1 Method

###  S1.1 Network details

We adopted a discrete-time formulation of network dynamics, in which

$\mathbf{x}_{t}=\mathbf{W}^{rec}\mathbf{r}_{t-1}+\mathbf{W}^{in}\mathbf{u}_{t}+\mathbf{W}^{in,att}[\mathbf{u}_{t}^{att}-\theta^{att}]_{+}+\mathbf{b}+\sqrt{2\sigma_{rec}^{2}}\text{N}(0,1),$
(S1)

where $\mathbf{x}_{t}$, $\mathbf{r}_{t}$ and $\mathbf{u}_{t}$ are respectively
the synaptic current, firing rate and network input at time step $t$,
$\mathbf{b}$ is the background input, $\mathbf{W}^{rec}$ is the recurrent
weight, $\mathbf{W}^{in}$ is the input weight, and $\sigma_{rec}$ is the
strength of recurrent noise. We supposed $\mathbf{r}_{t}=f(\mathbf{x}_{t})$,
with $f(\cdot)$ being the softplus current-rate transfer function, i. e.

$f(x)=\log(1+\exp(x)).$ (S2)

Input $\mathbf{u}_{t}$ is also noisy,

$\mathbf{u}_{t}=\mathbf{u}_{signal}+\sqrt{2\sigma_{in}^{2}}\text{N}(0,1),$
(S3)

with $\sigma_{in}$ being the strength of input noise. $\mathbf{W}^{in,att}$,
$\mathbf{u}_{t}^{att}$ and $\theta^{att}$ are the quantities related to the
input units modulated by top-down attention. They are only valid when studying
the effect of anticipatory attention in non-timing tasks (Fig. 6g-i). The
model does not have these quantities in the other tasks. $\mathbf{W}^{in,att}$
is the weight from the attention-modulated units to the recurrent network,
$\mathbf{u}_{t}^{att}$ is the input current to the attention-modulated units,
and $\theta^{att}$ is the firing threshold of these units. The firing
threshold is

$\theta^{att}=[\theta_{0}^{att}-\mathbf{W}^{fb,att}\mathbf{r}_{t}]_{+},$ (S4)

with $\mathbf{W}^{fb,att}$ being positive feedback weight, so that
$\theta^{att}$ decreases with feedback current until to zero, starting from
$\theta_{0}^{att}=1.5$. Eq. S4 models the disinhibitory effect of feedback
connections [1]. Similar to $\mathbf{u}_{t}$, $\mathbf{u}_{t}^{att}$ is also
noisy, with the noise strength $\sigma_{in}^{2}$ (eq. S3).

Some previous studies started with a continuous-time formulation, and obtained
the discrete-time version using Euler method (omitting the attention-modulated
units):

$\mathbf{x}_{t}=(1-\alpha)\mathbf{x}_{t-1}+\alpha(\mathbf{W}^{rec}\mathbf{r}_{t-1}+\mathbf{W}^{in}\mathbf{u}_{t}+\mathbf{b}+\sqrt{2\alpha^{-1}\sigma_{rec}^{2}}\text{N}(0,1)),$
(S5)

with $\alpha=\Delta t/\tau$ being the ratio of time step length $\Delta t$ and
membrane time constant $\tau$. In our study, we effectively set $\alpha=1$,
similarly as the scheme used in Ref. [2, 3]. We also set $\Delta t=20\text{
ms}$. The output of the network is supposed to be

$z=\mathbf{W}^{out}\mathbf{r}+\mathbf{b}^{out},$ (S6)

with the dimension of $z$ depending on tasks.

We set $\sigma_{in}=0.01$, $\sigma_{rec}=0.05$ when training the network.
After training, when plotting the neuronal activities in the perception epoch
(Fig. 2b-d), we kept $\sigma_{in}=0.01$, $\sigma_{rec}=0.05$ so that the
neuronal temporal profiles under different durations of perception epoch did
not fully overlap. When doing the other analysis, we turned off the noises by
default.

###  S1.2 Task details

####  S1.2.1 Timing tasks

Interval production task (IP). The network received from 2 input units: from
one came the two pulses that defined the time interval, and from the other
came the Go cue. The interval between the beginning of the simulation and the
onset of the first pulse was

$T_{start}\sim U(60\text{ ms},500\text{ ms}),$ (S7)

where $U(t_{1},t_{2})$ is a uniform distribution between $t_{1}$ and $t_{2}$.
The interval between the offset of the first pulse and the onset of the second
pulse was

$T\sim U(400\text{ ms},1400\text{ ms}).$ (S8)

Note that we set the range of $T$ to be $[400\text{ ms},1400\text{ ms}]$
during training, but after training, we only investigated the performance of
the network when $T\in[600\text{ ms},1200\text{ ms}]$. The reason is that
there were boundary effects if, after training, $T$ took a value close to 400
ms or 1400 ms: if $T$ was close to 400 ms, then the time interval produced by
the network was biased to be larger than $T$; whereas if $T$ was close to 1400
ms, then the produced interval was biased to be smaller than $T$. Such biases
were weak if $T$ took a middle value (Fig. S1e).

The interval between the offset of the second pulse and the onset of the Go
cue (i. e. , the delay period) was

$T_{delay}\sim U(600\text{ ms},1600\text{ ms}).$ (S9)

All input pulses (including the two pulses that defined the time interval, and
the Go cue) lasted for 60 ms, and had strength 1. Input units stayed at 0 when
there were no pulses.

The target output was a scalar. It stayed at zero from the beginning, jumped
to 1 at time $T$ after the offset of the Go cue, and kept at 1 until the end
of the simulation at 300ms afterwards.

Interval comparision task (IC). The network received two successive long-
lasting stimuli respectively from two input units. The first stimuli, which
came from the first unit, started at time $T_{start}$ after the beginning of
the simulation, and lasted for duration $T_{1}$. Then after a delay interval
$T_{delay}$ , the second stimuli, which came from the second unit, started,
and lasted for duration $T_{2}$.

$T_{start}\sim U(60\text{ ms},500\text{ ms}),\quad T_{1}\sim U(400\text{
ms},1400\text{ ms}),\quad T_{delay}\sim U(600\text{ ms},1600\text{ ms}),\quad
T_{2}\sim U(400\text{ ms},1400\text{ ms})$ (S10)

All the input stimuli had strength 1. Input units stayed at 0 when there were
no stimuli.

The target outputs were two scalars $\hat{z}_{0}$ and $\hat{z}_{1}$. Both
stayed at zero from the beginning. If $T_{1}>T_{2}$, then $\hat{z}_{0}$ jumped
to 1 at the offset of the second stimulus, and stayed at 1 until the end of
the simulation at 300 ms afterwards. Otherwise, $\hat{z}_{1}$ jumped to 1 at
the offset of the second stimulus.

Timed spatial reproduction task (t-SR). The network successively received
three pulses from three input channels. The first channel was a line that
coded spatial locations. This line contained 32 units, whose preferred
directions were uniformly spaced from -6 to 25. For unit $i$ with preferred
location $y_{i}$, its activity in a pulse with location $x$ was

$A_{in}(t)\exp[-\frac{1}{2}(\frac{|y_{i}-x|}{2})^{2}],$ (S11)

where $A_{in}(t)=1$ during the presentation of the pulse and $A_{in}(t)=0$ at
the other time. In our simulation, the spatial locations of the stimuli were
uniformly drawn from 0 to 19. The second and third channels were both scalar
inputs. The pulse from the second channel defined the time interval to be
remembered together with the pulse from the first channel. The pulse from the
third channel acted as Go cue. $T_{start}$, $T$ and $T_{delay}$ were
distributed similarly as in IP (eqs. S7-S9).

The target output was a line with 32 units, which represented response
location using similar tuning curves as the ones used for the input line (eq.
S11):

$\hat{z}_{i}=A_{out}(t)\exp[-\frac{1}{2}(\frac{|y_{i}-x|}{2})^{2}],$ (S12)

where the amplitude $A_{out}(t)$ stayed at zero from the beginning, jumped to
1 at time $T$ after the offset of the Go cue, and stayed at 1 until the end of
the simulation at 300 ms afterwards.

Timed decision making task (t-DM). The network received from three channels of
scalar inputs. From the first two channels came the stimuli whose strengths
were to be compared with each other, and from the last channel came the Go cue
pulse. Starting from the beginning of simulation, the first two channels were
set to 0 for duration $T_{start}$, and then jumped to $A_{1}$ and $A_{2}$
respectively; after $T$ time, these two channels were set to 0 again. The Go
cue pulse came at time $T_{delay}$ after the offset of the first two channels.
Here,

$A_{1}=\gamma+c,\quad A_{2}=\gamma-c,$ (S13)

where $\gamma$ was the average strength of these two stimuli and was
distributed as $\gamma\sim U(0.8,1.2)$, and $c$ measured the strength
difference of these two stimuli, and was distributed as

$c\sim U(\\{-0.08,-0.04,-0.02,-0.01,0.01,0.02,0.04,0.08\\}),$ (S14)

where $U(\\{a_{1},a_{2},\cdots,a_{n}\\})$ denotes a discrete uniform
distribution over the set $\\{a_{1},a_{2},\cdots,a_{n}\\}$. $T_{start}$, $T$
and $T_{delay}$ were distributed similarly as in interval production task
(eqs. S7-S9).

The target outputs were two scalars $\hat{z}_{0}$ and $\hat{z}_{1}$. Both
stayed at zero from the beginning. If $c>0$, then $\hat{z}_{0}$ jumped to 1 at
time $T$ after the offset of the Go cue, and stayed at 1 until the end of the
simulation at 300ms afterwards. Otherwise, $\hat{z}_{1}$ jumped to 1 at time
$T$ after the offset of the Go cue.

####  S1.2.2 Non-timing tasks: default settings

Spatial reproduction task (SR). The network received pulses from two input
channels. The first channel was a line that contained 32 units, coding spatial
locations in the range $[-6,25]$ in the way indicated by eq. S11. In our
simulation, the spatial locations of the stimuli were uniformly drawn from 0
to 19. The second channel is a scalar input. The duration $T_{delay}$ of the
delay epoch between the first and second pulses was 1200 ms. The target output
was a line of 32 units (eq. S12), which was to indicate the location of the
first pulse immediately after the second pulse.

Comparison task (COMP). The network received pulses from two input channels,
both of which were lines that contained 32 units successively gave two pulses
to the network. The target outputs were two scalars $\hat{z}_{0}$ and
$\hat{z}_{1}$, which were to indicate whether or not the spatial coordinate of
the first pulse was larger than that of the second pulse.

Change detection task (CD). The network had the same structure as that in
COMP. Two scalar outputs were to indicate whether or not the distance between
the spatial locations of the two input pulses was within 1.

Decision making task (DM). The network received from two channels of stimuli
lasting for $T=1200$ ms. The two scalar outputs were to indicate which
stimulus was stronger immediately after the ending of the two stimuli.

Cue-dependent decision making task (cue-DM). The network received from two
channels of stimuli lasting for $T=1200$ ms. At the 1140 ms after the
presentation of the two stimuli, a two dimensional one-hot vector lasting for
60ms was input from a third channel. Two scalar outputs were to indicate the
index of the stronger stimulus or the index of the weaker stimulus according
to the third channel.

####  S1.2.3 Non-timing tasks: studying the factors that influence the
strength of temporal signal

To study the effect of the overlap of sensory input to the strength of
temporal signal in the delay epoch of SR, COMP and CD (Fig. 6d), we expanded
the unit number in the line channels to 44 (default is 32), and broadened the
standard deviation of the tuning curves (eq. S12) to 4 (default is 2). These
units coded spatial locations in the range -12 to 31. In our simulation, the
spatial locations of input stimuli were uniformly drawn from 0 to 19.

To study the effect of multi-tasking (Fig. 6e), we trained the network on t-SR
and SR concurrently, or on t-DM and DM concurrently. The two tasks in each
pair share the same input and output channels. We used a one-hot vector from
another two-dimensional input channel to indicate which task should be
performed [4]. The network was to be able to perform either of the indicated
task.

To study the effect of timing anticipation (Fig. 6f), we trained the network
to perform SR, COMP, CD, DM and cue-DM, with the duration $T$ of the delay
epoch (for SR, COMP and CD) or the stimuli-presentation epoch (for DM and cue-
DM) was randomly between $[800\text{ ms},1600\text{ ms}]$. After training, we
analyzed the simulation results when $T=1200$ ms, and compared the results
with the cases after training the network with $T$ fixed at 1200 ms. To study
the the effect of anticipatory attention (Fig. 6g-i), feedback was imposed on
the second input channel of SR, COMP and CD, and was imposed on the third
channel of cue-DM. This means that these input channels were modeled using the
third term at the right-hand side in eq. S1, instead of the second term.

###  S1.3 Training details

Training was performed to minimize a cost function using back-propagation
through time. Cost function was defined as

$C=\sum_{i}m_{i}(z_{i}-\hat{z}_{i})^{2},$ (S15)

where $i$ is the index of output units, $z_{i}$ is the actual output defined
by eq. S6, $\hat{z}_{i}$ is the target output, and $m_{i}$ is the mask. In all
tasks, $m_{i}=0$ before the onset of the first stimulus, and $m_{i}=1$
afterwards; therefore, only the output after the onset of the first stimulus
was constrained. When studying the effect of anticipatory attention in non-
timing tasks (Fig. 6g-i), we added L2 regularization to feedback current
$\mathbf{I}^{fb}=\mathbf{W}^{fb,att}\mathbf{r}_{t}$ (see eq. S4), so that
eq.S15 becomes
$C=\sum_{i}m_{i}(z_{i}-\hat{z}_{i})^{2}+\beta_{fb}\frac{1}{N_{i,t}}\sum_{i,t}(I_{i,t}^{fb})^{2}$,
with $\beta_{fb}=10^{-4}$. This cost function was minimized using Adam
optimizer at learning rate 0.0005, with batch size 64 in each training step.
We trained 16 configurations to perform IP and IC tasks, and trained 30
configurations to perform t-SR and t-DM tasks. Different configurations were
initialized using different random seeds.

Before training, recurrent self-connections ($W_{ii}^{rec}$ in eq. S5) were
initialized to 1, and other recurrent connections were initialized as
independent Gaussian variables with mean 0 and standard deviation
$0.3/\sqrt{N_{rec}}$, with $N_{rec}=256$ being the number of recurrent units.
This initialization strategy was used in Ref. [3]. The identity self-
connections prevent vanishing gradient during training [5], and the non-zero
off-diagonal recurrent connections induce sequential activity in the network
after training [3], so that the dynamics of the network becomes comparable to
experimental observations [6, 7, 8, 9, 10]. Output connections were
initialized as independent Gaussian variable with mean 0 and standard
deviation $1/\sqrt{N_{rec}}$. Input connections from the line input were
initialized as variables drawn uniformly from
$[-1/\sqrt{2\sigma_{tuning}},1/\sqrt{2\sigma_{tuning}}]$, with
$\sigma_{tuning}$ being the standard deviation of the Gaussian tuning curve
(eq. S11), which was 2 by default and 4 when studying the effect of input
overlap in non-timing tasks. The input connections from the other channels
were initialized as variables drawn uniformly from
$[-1/\sqrt{D_{channel}},1/\sqrt{D_{channel}}]$, with $D_{channel}$ being the
dimension of the input channel.

Every 200 training steps, we evaluated the performance of the network using a
batch of size 512, and stopped training as soon as the performance of the
network reached criterion (Fig. S1i-l). We introduced our criterion in t-SR
and t-DM in details, the other tasks shared similar criterion:

In t-SR, a time interval was considered to be produced if: (1) the activities
of all the 32 output units were below 0.2 before the offset of the Go cue, (2)
one of them went above 0.5 at some point $t_{p}$ before $T+300\text{ms}$ after
the offset time $t_{off}^{cue}$ of the Go cue. The produced interval was
$T_{p}=t_{off}^{cue}-t_{p}$. Output location at time $t_{p}$ was read out
using a population vector method (see the computer code in Ref. [4]). Training
was stopped as soon as (1) time intervals were produced in over 95% simulation
trials, (2) the relative error of the produced intervals $|T_{p}-T|/T<0.025$,
(3) the output locations were on average within 0.8 of the input locations.

In t-DM, a time interval was considered to be produced if: (1) the activities
of both output units $z_{0}$ and $z_{1}$ were below 0.2 before the offset of
the Go cue, (2) one of them went above 0.5 at some time point $t_{p}$ before
$T+300\text{ms}$ after the offset $t_{off}^{cue}$ of the Go cue, whereas the
other one stayed below 0.5. The produced interval was
$T_{p}=t_{off}^{cue}-t_{p}$. In the trials in which a time interval was
produced, the decision was considered to be correct if: when $c>0$ (or $c<0$),
$z_{0}$ (or $z_{1}$) went above 0.5 and $z_{1}$ (or $z_{0}$) kept below 0.5.
Training was stopped as soon as (1) time intervals were produced in over 96%
of simulation trials, (2) the relative error of the produced intervals
$|T_{p}-T|/T<0.025$, (3) the decision error rate was smaller than 0.02.

###  S1.4 Data analysis

####  S1.4.1 Types of neurons at the end of the delay epoch

In IP or IC, we supposed $f_{i}(T)$ to be the activity of the $i$th neuron at
the end of the delay epoch as a function of the duration $T$ of the perception
(for IP) or stimulus1 (for IC) epoch. We picked neurons that can be strongly
activated at the end of the delay epoch, namely the neurons whose
$\max_{T\in[T_{min},T_{max}]}f_{i}(T)>\theta_{sa}$, with $T_{min}=600\text{
ms}$ and $T_{max}=1200\text{ ms}$ respectively being the minimal and maximal
values of $T$ in our simulation, and $\theta_{sa}=2$. Our results are not
sensitive to the value of $\theta_{sa}$. We classified $f_{i}(T)$ of the
picked neurons into three types, namely monotonically increasing (MoI),
monotonically decreasing (MoD), and non-monotonic (non-M) in the following
way: We divided the range of $T$ (i.e., $[T_{min},T_{max}]$) into four parts
of the same length, and calculated the mean value of $f_{i}(T)$ in these four
parts, say $f_{i}(\text{part
1})=\frac{4}{T_{max}-T_{min}}\int_{T_{min}}^{T_{min}+(T_{max}-T_{min})/4}f_{i}(T)\mathrm{d}T$,
$f_{i}(\text{part
2})=\frac{4}{T_{max}-T_{min}}\int_{T_{min+(T_{max}-T_{min})/4}}^{T_{min}+2(T_{max}-T_{min})/4}f_{i}(T)\mathrm{d}T$,
etc. If $f_{i}(\text{part 1})\leq f_{i}(\text{part 2})\leq f_{i}(\text{part
3})\leq f_{i}(\text{part 4})$, then neuron $i$ belongs to MoI type; if
$f_{i}(\text{part 1})\geq f_{i}(\text{part 2})\geq f_{i}(\text{part 3})\geq
f_{i}(\text{part 4})$, then neuron $i$ belongs to MoD type; otherwise, neuron
$i$ belongs to non-M type.

In t-SR, we supposed $g_{i}(T,x)$ to be the activity of the $i$th neuron at
the end of the delay epoch as a function of $T$ at a given location $x$ of the
first pulse. We picked neurons that can be strongly activated at the end of
the delay epoch (i.e., the neurons whose
$\max_{\\{T,x\\}}g_{i}(T,x)>\theta_{sa}$). We then defined
$f_{i}(T)=\max_{x}g_{i}(T,x)$, and classified neuron $i$ into MoI, MoD or
non-M types according to the monotonicity of $f_{i}(T)$ in the similar way to
the IP or IC case introduced above. Similarly, in t-DM, we classified neurons
according to $f_{i}(T)=\max_{c}g_{i}(T,c)$, where $c$ is the half difference
between the strengths of the presented stimuli (eq. S13).

####  S1.4.2 Temporal scaling in the production epoch

Analysis of temporal scaling was performed using similar technique to Ref.
[2]. Specifically, we calculated the $k$th scaling component
$\mathbf{u}_{SC,k}$ through the following equation:

$\mathbf{u}_{SC,k}=\text{arg
}\min_{\mathbf{u}}\frac{\sum_{t}\sum_{T}(\mathbf{r}_{k}^{S}(t;T)\mathbf{u}-\text{Mean}_{T}(\mathbf{r}_{k}^{S}(t;T)\mathbf{u}))^{2}}{\sum_{t}\sum_{T}(\mathbf{r}_{k}^{S}(t;T)\mathbf{u}-\text{Mean}_{\\{t,T\\}}(\mathbf{r}_{k}^{S}(t;T)\mathbf{u}))^{2}},$
(S16)

where $\mathbf{r}_{k}^{S}(t;T)$ is population activity at the scaled time when
the duration of the perception epoch is $T$ (see below for details), the
denominator is the total variance of the trajectories, and the numerator is
the variance that cannot be explained by temporal scaling. To calculate the
first scaling component $\mathbf{u}_{SC,1}$, we set
$\mathbf{r}_{1}^{S}(t;T)=\mathbf{r}^{PC}(tT_{p};T),$ with $0\leq t\leq 1$,
where $\mathbf{r}^{PC}$ is the projection of the population activity in the
subspace spanned by the first 9 principal components, and $T_{p}$ is the
interval produced by the network in the production epoch; then we minimized
$\mathbf{u}$ in eq. S16. To calculate the second scaling component
$\mathbf{u}_{SC,2}$, we set
$\mathbf{r}_{2}^{S}(t;T)=\mathbf{r}_{1}^{S}(t;T)-\mathbf{r}_{1}^{S}(t;T)\mathbf{u}_{SC,1}$,
and then minimized $\mathbf{u}$ in eq. S16 in the subspace orthogonal to
$\mathbf{u}_{SC,1}$. In this way, we calculated all the 9 scaling components
one by one.

Scaling index (SI) of a subspace $U$ was defined as

$\text{SI}=\frac{\sum_{t}\sum_{T}(\mathbf{r}_{1}^{S}(t;T)U-\text{Mean}_{T}(\mathbf{r}_{1}^{S}(t;T)U))^{2}}{\sum_{t}\sum_{T}(\mathbf{r}_{1}^{S}(t;T)U-\text{Mean}_{\\{t,T\\}}(\mathbf{r}_{1}^{S}(t;T)U))^{2}},$
(S17)

where $\mathbf{r}_{1}^{S}(t;T)U$ is the projection of the scaled trajectory to
the subspace $U$.

####  S1.4.3 The geometry of coding combination

During the perception epoch of t-SR, the network state is quantified by the
time elapsed from the beginning of the epoch (temporal flow) and the spatial
information of the first pulse. At the end of the delay epoch of t-SR, the
network state is quantified by the time interval between the first two pulses
and the spatial information of the first pulse. During the production epoch of
t-SR, the network state is quantified by temporal flow, time interval and
spatial information. Similar scenario also exists in t-DM, except that the
non-temporal information is the decision choice made by the network. In t-DM,
the decision choice $d$ depends on the sign of the half difference $c$ between
the strength of the presented two stimuli (eq. S13), we defined
$r_{i}(d=1,\\{a\\})=\langle r_{i}(c,\\{a\\})\rangle_{c>0}$ and
$r_{i}(d=-1,\\{a\\})=\langle r_{i}(c,\\{a\\})\rangle_{c<0}$, where $\\{a\\}$
indicates the other parameters than decision choice, and used
$r_{i}(d,\\{a\\})$ to do the following analysis. Together, during the
perception epoch and at the end of the delay epoch of t-SR and t-DM, two
variables are coded in the network state; during the production epoch, three
variables are coded in the network state. We used two measurements to quantify
the geometry of the coding combination of multiple variables: (1) the angle
between the first marginal principal components and (2) the mixed variance
[11], introduced below.

Suppose the activity of the $i$th neuron $r_{i}(a,b)$ is a function of two
variables $a$ and $b$, with the mean of $r_{i}(a,b)$ being subtracted so that
$\langle r_{i}(a,b)\rangle_{a,b}=0$. The marginal principal components (PCs)
with respect to $a$ are the PCs of the dot set $\\{\langle
r_{i}(a,b)\rangle_{b}\\}_{i}$, and the marginal PCs of $b$ are the PCs of
$\\{\langle r_{i}(a,b)\rangle_{a}\\}_{i}$. We quantified the coding
orthogonality of $a$ and $b$ by calculating the angle between the first
marginal PCs of $a$ and $b$. The portions of variance explained by $a$ and $b$
are respectively $p_{a}=\text{Var}_{i,a}(\\{\langle
r_{i}(a,b)\rangle_{b}\\}_{i})/v_{tot}$ and $p_{b}=\text{Var}_{i,b}(\\{\langle
r_{i}(a,b)\rangle_{a}\\}_{i})/v_{tot}$, with the total variance
$v_{tot}=\text{Var}_{i,a,b}(\\{r_{i}(a,b)\\}_{i})$. The portion of mixed
variance between $a$ and $b$ is $p_{a+b}=1-p_{a}-p_{b}$.

In the case that the activity of the $i$th neuron $r_{i}(a,b,c)$ is a function
of three variables, we also subtracted the mean of $r_{i}(a,b,c)$ so that
$\langle r_{i}(a,b,c)\rangle_{a,b,c}=0$. The marginal PCs of $a$, $b$ and $c$
are respectively the PCs of $\\{\langle r_{i}(a,b,c)\rangle_{b,c}\\}_{i}$,
$\\{\langle r_{i}(a,b,c)\rangle_{a,c}\\}_{i}$ and $\\{\langle
r_{i}(a,b,c)\rangle_{a,b}\\}_{i}$. The portions of variance explained by these
variables and their mixing were defined as [11]:

$p_{a}=\text{Var}_{i,a}(\\{\langle r_{i}(a,b)\rangle_{b,c}\\}_{i})/v_{tot}$
$p_{b}=\text{Var}_{i,b}(\\{\langle r_{i}(a,b)\rangle_{a,c}\\}_{i})/v_{tot}$
$p_{c}=\text{Var}_{i,c}(\\{\langle r_{i}(a,b)\rangle_{a,b}\\}_{i})/v_{tot}$
$p_{a+b}=\text{Var}_{i,a,b}(\\{\langle r_{i}(a,b,c)-\langle
r_{i}(a,b)\rangle_{b,c}-\langle r_{i}(a,b)\rangle_{a,c}-\langle
r_{i}(a,b)\rangle_{a,b}\rangle_{c}\\}_{i})/v_{tot}$
$p_{b+c}=\text{Var}_{i,b,c}(\\{\langle r_{i}(a,b,c)-\langle
r_{i}(a,b)\rangle_{b,c}-\langle r_{i}(a,b)\rangle_{a,c}-\langle
r_{i}(a,b)\rangle_{a,b}\rangle_{a}\\}_{i})/v_{tot}$
$p_{a+c}=\text{Var}_{i,a,c}(\\{\langle r_{i}(a,b,c)-\langle
r_{i}(a,b)\rangle_{b,c}-\langle r_{i}(a,b)\rangle_{a,c}-\langle
r_{i}(a,b)\rangle_{a,b}\rangle_{b}\\}_{i})/v_{tot}$
$p_{a+b+c}=1-p_{a}-p_{b}-p_{c}-p_{a+b}-p_{b+c}-p_{a+c}$

where $v_{tot}=\text{Var}_{i,a,b,c}(\\{r_{i}(a,b,c)\\}_{i})$ is the total
variance, “$+$” sign in the subscript indicates the mixing of several
variables.

In Fig. 3, we used the network state trajectory after 400 ms (200 ms) of
transient period of the perception (production) epoch to do the analysis.

####  S1.4.4 Decoding

We studied two types of nearest-centroid decoders [12]. Given a population
state $\mathbf{f}_{0}$, the decoded value $a_{d,1}$ read-out by Decoder 1 is

$a_{d,1}=\text{arg
min}_{a\in\mathcal{A}}(\left\|\mathbf{f}_{0}\mathbf{W}^{dec}-\mathbf{f}(a;b_{train})\mathbf{W}^{dec}\right\|),$
(S18)

where $\mathbf{f}(a;b_{train})$ is the population state as a function of
variable $a$ along an iso-$b$ line whose $b$ value is constantly $b_{train}$,
and decoding weight $\mathbf{W}^{dec}$ is the first PC of
$\mathbf{f}(a;b_{train})$. The decoded value $a_{d,2}$ read-out by Decoder 2
is

$a_{d,2}=\text{arg
min}_{a\in\mathcal{A}}(\left\|(\mathbf{f}_{0}-\langle\mathbf{f}(a;b_{test})\rangle_{a})\mathbf{W}^{dec}-(\mathbf{f}(a;b_{train})-\langle\mathbf{f}(a;b_{train})\rangle_{a})\mathbf{W}^{dec}\right\|),$
(S19)

where $\mathbf{f}(a;b_{test})$ is the iso-$b$ line that $\mathbf{f}_{0}$
belongs to, and $\langle\cdot\rangle_{a}$ means averaging over $a$. From
eq.S19, both the mass centers of the two iso-$b$ lines
$\mathbf{f}(a;b_{train})$ and $\mathbf{f}(a;b_{test})$ are translationally
moved to the zero point before $\mathbf{f}(a;b_{train})$ and
$\mathbf{f}(a;b_{test})$ are projected to the decoding space by
$\mathbf{W}^{dec}$.

####  S1.4.5 Correlation between decoding error, angle and mixed variance

In Fig. 4d, f, we computed the correlation between decoding error (DE), the
angle (AG) between the first PCs of the decoded and generalized variables, and
the mixed variance (MV) between the decoded and generalized variables. A
subtle point here is that AG and MV may also be correlated (see Fig. 4c, e for
the negative correlation between AG and MV in the production epoch of t-SR),
therefore the Pearson’s correlation between DE and AG may be contributed by
two pathways: (1) AG influences DE directly; (2) AG influences DE indirectly
through MV, due to the correlation between AG and MV. Similar situation also
exists for the correlation between DE and MV. To investigate the direct
correlation and remove the indirect one, we iteratively took the following
operation to reduce the correlation between AG and MV: removing a single data
point (i.e., the AG and MV of a single training configuration) from the
dataset, so that the absolute value of the correlation between AG and MV in
the left dataset is minimal. We found that small correlation (with absolute
value below 0.05) between AG and MV could usually be obtained after removing 2
or 3 data points from the whole dataset of 30 points (Figs. S7, S8). In this
way, we got a dataset with small correlation between AG and MV, while at the
same time, as large as possible. Pearson’s correlation were then calculated
using the left dataset to draw Figs. 4d, f, S7, S8.

####  S1.4.6 Firing sequence and network structure

To plot Fig. 5a, b, we ordered the peak firing time of strongly active neurons
(whose peak firing rates were larger than 2) in the studied epoch, and plotted
weight connection as a function of the peak order difference between the post-
and pre-synaptic neurons.

To plot Fig. 5c, d, we used a more elaborate method to illustrate the network
structure underlying t-SR and t-DM. At time $t_{0}$ and non-time information
$x_{0}$ (which may be spatial location or decision choice), we picked a set
$\mathcal{N}(t_{0},x_{0})$ of strongly active neurons whose firing rates at
$t_{0}$ and $x_{0}$ were larger than a threshold 2 (our result is insensitive
to this threshold). We then defined $T_{peak,i}(t_{0},x_{0})$ to be the peak
time of neuron $i$ near $t_{0}$ at $x_{0}$: if the activity $f_{i}(t;x_{0})$
of neuron $i$ decreased (or increased) with time at time point $t_{0}$ and
non-time information $x_{0}$, then $T_{peak,i}(t_{0},x_{0})$ was the time
point of the local maximum of $f_{i}(t;x_{0})$ before (or after), but most
nearest to, $t_{0}$. Iterating over all the possible values of $x_{0}$, we got
all the strongly active neurons at time $t_{0}$:
$\mathcal{N}(t_{0})=\bigcup_{x_{0}}\mathcal{N}(t_{0},x_{0})$. For neuron $i$
in $\mathcal{N}(t_{0})$, we called its prefered non-time information
$x_{prefer}$ to be the value of $x_{0}$ that maximized its peak firing rate:
$x_{prefer}=\text{arg max}_{x_{0}}f_{i}(T_{peak,i}(t_{0},x_{0}),x_{0})$. In
this way, we classified all the neurons in $\mathcal{N}(t_{0})$ according to
their non-time information preference:
$\mathcal{N}(t_{0})=\bigcup_{x_{0}}\mathcal{N}_{prefer}(t_{0},x_{0})$, with
$\mathcal{N}_{prefer}(t_{0},x_{0})$ being the set of neurons that prefer
$x_{0}$ around time $t_{0}$. We then defined $T_{peak,i}(t_{0},x_{prefer})$ to
be the big peak time of neuron $i$ at time $t_{0}$. Given a neuron $i$ and a
set $\mathcal{N}_{prefer}(t_{0},x_{0})$ of neurons ($i$ may or may not belong
to $\mathcal{N}_{prefer}(t_{0},x_{0})$), we ordered their big peak times, and
then investigated the recurrent weight from $i$ to each neuron of
$\mathcal{N}_{prefer}(t_{0},x_{0})$ (except $i$ itself if
$i\in\mathcal{N}_{prefer}(t_{0},x_{0})$). In this way, we studied the
recurrent weight $w(o_{post}-o_{pre},|x_{post}-x_{pre}|)$ as a function of the
difference $o_{post}-o_{pre}$ between the orders of the big peak time of the
post- and pre-synaptic neurons and the difference $|x_{post}-x_{pre}|$ of
their preferred non-time information. Fig. 5c,d were plotted by averaging
$w(o_{post}-o_{pre},|x_{post}-x_{pre}|)$ over $t_{0}$ and training
configurations.

##  S2 The relationship between the low dimensionality of the attractor in the
delay epoch and the dominance of monotonic neurons

We denote $\mathcal{M}$ as the manifold of the population states at the end of
the delay epoch at different durations $T$ of the perception epoch (Fig. 2e).
The first principal component (PC) of $\mathcal{M}$ explained about 90% of its
variance (Fig. 2g), and the activities of most neurons changed monotonically
with $T$ in $\mathcal{M}$ (Fig. 2j). To understand the relationship between
these two facts, let’s consider the extreme case that all neurons are linearly
monotonic with $T$ in $\mathcal{M}$, then $\mathcal{M}$ is a line in the
population-state space that can be parameterized as
$[f_{1}(T),f_{2}(T),\cdots,f_{N}(T)]^{T}$, with $f_{i}(T)$ being the activity
of the $i$th neuron at the end of the delay epoch when the duration of the
perception epoch is $T$. In this case, PC1 of $\mathcal{M}$, which explains
100% of the variance of $\mathcal{M}$ because $\mathcal{M}$ is a line, is the
following vector with unit length:

$\pm\frac{1}{\sqrt{\sum_{i}(f_{i}(T_{max})-f_{i}(T_{min}))^{2}}}[f_{1}(T_{max})-f_{1}(T_{min}),f_{2}(T_{max})-f_{2}(T_{min}),\cdots,f_{N}(T_{max})-f_{N}(T_{min})]^{T},$

where $T_{min}=600\text{ms}$ and $T_{max}=1200\text{ms}$ are respectively the
minimal and maximal values of $T$ in our simulation, and the $\pm$ sign
indicates that the direction of PC1 is undetermined. If neuron $i$
monotonically increases (or decreases) with $T$, then
$f_{i}(T_{max})-f_{i}(T_{min})>0$ (or $f_{i}(T_{max})-f_{i}(T_{min})<0$).
Apparently, if two neurons $i$ and $j$ have the same (or different)
monotonicity, then their corresponding elements in PC1 have the same
(different) signs. This is indeed what we found in our simulation (Fig. S2g,
h).

##  S3 The geometric meaning of mixed variance

We denote the population state to be
$\mathbf{r}=\\{r_{1},r_{2},\cdots,r_{N}\\}$, where $r_{i}$ is the firing rate
of the $i$th neuron, or in general, the activity projected on the $i$th basis
vector, say, principal component. Suppose $\mathbf{r}$ is parameterized by two
variables $a$ and $b$, and we subtract the mean value of $r_{i}$ so that

$\text{E}_{a,b}[r_{i}(a,b)]=0,$ (S20)

where $\text{E}_{a,b}[\cdot]$ means the average over $a$ and $b$.

The total variance of $\mathbf{r}$ is

$v_{tot}=\text{Var}_{i,a,b}[r_{i}(a,b)]$
$=\text{E}_{i}[\text{Var}_{a,b}[r_{i}(a,b)]]+\text{Var}_{i}[\text{E}_{a,b}[r_{i}(a,b)]]$
$=\text{E}_{i}[\text{Var}_{a,b}[r_{i}(a,b)]],$ (S21)

where $\text{Var}_{x}[\cdot]$ means the variance over variable $x$. The first
equation is the definition of the total variance, the second equation is from
the law of total variance, and the third equation is from eq. S20. Similarly,
the variance explained by $a$ is

$v_{a}=\text{Var}_{i,a}[\text{E}_{b}[r_{i}(a,b)]]=\text{E}_{i}[\text{Var}_{a}[\text{E}_{b}[r_{i}(a,b)]]],$
(S22)

and the variance explained by $b$ is

$v_{b}=\text{Var}_{i,b}[\text{E}_{a}[r_{i}(a,b)]]=\text{E}_{i}[\text{Var}_{b}[\text{E}_{a}[r_{i}(a,b)]]]$
(S23)

Now let’s study a sufficient condition so that

$v_{tot}=v_{a}+v_{b},$ (S24)

which means that the mixed variance

$v_{mix}=v_{tot}-(v_{a}+v_{b})$ (S25)

is zero.

From eqs. S21-S23, a sufficient condition to fulfill eq. S24 is

$\text{Var}_{a,b}[r_{i}(a,b)]=\text{Var}_{a}[\text{E}_{b}[r_{i}(a,b)]]+\text{Var}_{b}[\text{E}_{a}[r_{i}(a,b)]]\quad\text{for
every }i.$ (S26)

According to the law of total variance,

$\text{Var}_{a,b}[r_{i}(a,b)]=\text{Var}_{a}[\text{E}_{b}[r_{i}(a,b)]]+\text{E}_{a}[\text{Var}_{b}[r_{i}(a,b)]].$
(S27)

Therefore, to realize eq. S26, we can set

$\text{Var}_{b}[\text{E}_{a}[r_{i}(a,b)]]=\text{E}_{a}[\text{Var}_{b}[r_{i}(a,b)]]\quad\text{for
every }i.$ (S28)

in other words

$\text{E}_{b}[(\text{E}_{a}[r_{i}(a,b)]-\text{E}_{a,b}[r_{i}(a,b)])^{2}]=\text{E}_{a}[\text{E}_{b}[(r_{i}(a,b)-\text{E}_{b}[r_{i}(a,b)])^{2}]]\quad\text{for
every }i$ (S29)

Because $\text{E}_{a,b}[r_{i}(a,b)]=0$, this equation gives

$\text{E}_{b}[(\text{E}_{a}[r_{i}(a,b)])^{2}]=\text{E}_{a}[\text{E}_{b}[(r_{i}(a,b)-\text{E}_{b}[r_{i}(a,b)])^{2}]]\quad\text{for
every }i$ (S30)

A sufficient condition to fulfill the equation above is

$r_{i}(a,b)-\text{E}_{b}[r_{i}(a,b)]=f(b)\quad\text{for every }i,$ (S31)

namely the value of $r_{i}(a,b)-\text{E}_{b}[r_{i}(a,b)]$ does not depend on
$a$. This sufficient condition can be easily proved by substituting eq. S31
into eq. S30 and using the fact that $\text{E}_{a,b}[r_{i}(a,b)]=0$. Now let’s
try to understand the meaning of eq. S31. Consider four pairs of variables
$(a_{1},b_{1})$, $(a_{2},b_{1})$, $(a_{1},b_{2})$ and $(a_{2},b_{2})$, we have

$r_{i}(a_{1},b_{1})-\text{E}_{b}[r_{i}(a_{1},b)]=f(b_{1})=r_{i}(a_{2},b_{1})-\text{E}_{b}[r_{i}(a_{2},b)]\quad\text{for
every }i$ (S32)
$r_{i}(a_{1},b_{2})-\text{E}_{b}[r_{i}(a_{1},b)]=f(b_{2})=r_{i}(a_{2},b_{2})-\text{E}_{b}[r_{i}(a_{2},b)]\quad\text{for
every }i$ (S33)

By subtracting eq. S32 from eq. S33, we have

$r_{i}(a_{1},b_{1})-r_{i}(a_{1},b_{2})=r_{i}(a_{2},b_{1})-r_{i}(a_{2},b_{2})\quad\text{for
every }i.$ (S34)

This means that between the two iso-$b$ lines in which the values of $b$ are
separately fixed at $b_{1}$ and $b_{2}$, the vector that connects the two
points representing $a_{1}$ is equal to the vector that connects the two
points representing $a_{2}$. In other words, these two iso-$b$ lines can be
related by translational movement. By rewritten eq.S34 as
$r_{i}(a_{1},b_{1})-r_{i}(a_{2},b_{1})=r_{i}(a_{1},b_{2})-r_{i}(a_{2},b_{2})$,
we see that different iso-$a$ lines are also related by translational
movement.

From the discussion above, translational relation between different iso-$a$ or
iso-$b$ lines is a sufficient condition for zero mixed variance. How about the
necessity? In other words, if we observe close-to-zero mixed variance in
simulation, how will be the geometry of the iso-$a$ and iso-$b$ lines? We
checked this point through simulation. In Fig. S6, we show the iso-space lines
of several simulation examples, in the perception, delay and production epochs
of t-SR task. We see that in examples with small mixed variance, the iso-space
lines of different spatial information tend to be parallel and of the same
length; whereas in examples with large mixed variance, the iso-space lines may
be non-parallel or of very different lengths. Additionally, if iso-$a$ or
iso-$b$ lines are translationally related, then Decoder 2 (eq. S19) will have
perfectly zero generalization error. We found that the generalization error of
Decoder 2 is strongly positively correlated with mixed variance (Figs. 4f, S7,
S8). These results imply that at least in the context of our simulation, mixed
variance is a good index to quantify the translational relationship between
different iso-$a$ or iso-$b$ lines, or in other words, the parallelogram-
likeness of iso-$a$ and iso-$b$ grids (Fig. 3f, upper left).

The opposite extreme case that $v_{mix}=v_{tot}$, which, from eq.S25, means
$v_{a}=v_{b}=0$. From eqs. S22, S23, this means that

$\text{Var}_{a}[\text{E}_{b}[r_{i}(a,b)]]=\text{Var}_{b}[\text{E}_{a}[r_{i}(a,b)]]=0\quad\text{for
every }i.$

In other words, the mean value of $r_{i}(a,b)$ over $b$ (i.e.,
$\text{E}_{b}[r_{i}(a,b)]$) does not depends on $a$, and the mean value of
$r_{i}(a,b)$ over $a$ (i.e., $\text{E}_{a}[r_{i}(a,b)]$) does not depends on
$b$ neither. This implies that different iso-$a$ (and also iso-$b$) lines are
strongly intertwined with each other, so that they have the same mean state
value. A good example of this case is that every point in the 2-dimensional
range of variables $[a_{min},a_{max}]\otimes[b_{min},b_{max}]$ (where
$a_{min}$, $a_{max}$, $b_{min}$ and $b_{max}$ are the minimal and maximal
values of $a$ and $b$ respectively) is mapped toward a random point in a state
space
$[r_{1,min},r_{1,max}]\otimes[r_{2,min},r_{2,max}]\otimes\cdots\otimes[r_{n,min},r_{n,max}]$:
in this case, every iso-$a$ or iso-$b$ dot set of states has the mean value
located at the center of the state space
$(\frac{r_{1,min}+r_{1,max}}{2},\frac{r_{2min}+r_{2,max}}{2},\cdots,\frac{r_{n,min}+r_{n,max}}{2})$.

Figure S1: Performance of the network after training. (a-e) Interval
production (IP) task. (a) An example of the input and output of the network in
IP. Red and blue lines: two input channels. Dashed black line: target output.
Solid black line: actual output. (b) Probability distribution function (p.d.f)
of self-connections (blue) and non-diagonal connections (red) of the recurrent
network after training. (c) Three examples of the output in the production
epoch of IP, when $T=600$ ms (blue), 900 ms (red) and 1200 ms (yellow). Dashed
line: target output. Solid line: actual output. The horizontal dashed black
line indicates the threshold that the network is regarded to generate a
movement in the production epoch when the output rises across this threshold.
(d) Distribution of the scaling index of the output across training
configurations in the production epoch of IP. (e) The difference between the
produced time interval $T_{p}$ and the interval $T$ between the first two
pulses in IP as a function of $T$. Error bar means standard deviation over 16
training configurations. During training, we set $T\in[400\text{
ms},1400\text{ ms}]$. This panel shows that if after training we set $T$ to be
close to 400 ms, $T_{p}$ tends to be larger than $T$; whereas if we set $T$ to
be close to 1400 ms, $T_{p}$ tends to be smaller than $T$. Therefore, by
default, we set $T\in[600\text{ ms},1200\text{ ms}]$ for data analysis after
training to reduce the bias of $T_{p}$. (f) Two examples of interval
discrimination (IC) task. Upper: the case when the duration of the first
stimulus is shorter than that of the second stimulus. Lower: the case when the
duration of the first stimulus is longer than that of the second stimulus. Red
and yellow lines: two input channels. Dashed black and pink lines: two
channels of target output. Solid black and pink lines: two channels of actual
output. (g) An example of timed spatial reproduction (t-SR) task. Left upper:
the pulse with location information from the first input channel. Left lower:
the pulses from the second (yellow) and third (blue) input channels. Right:
actual output. (h) Two examples of timed decision making (t-DM) task. Upper:
when the input from the first channel (red) is weaker than the input from the
second channel (yellow), i.e., $c<0$. Lower: when $c>0$. (i-l) Performance of
the network during training. (i) Performance of the network during the
training of IP, quantified by the probability to successfully produce time
interval (upper) and the relative error of the produced interval (lower). Gray
lines indicate individual training configurations. Training stopped as soon as
both quantities reach the criterion (horizontal dashed lines). (j) Performance
of the network during the training of IC, quantified by the probability to
successfully output a choice (upper) and the probability of choice error
(lower). (k) Performance of the network during the training of t-SR,
quantified by the probability to successfully produce time interval (upper),
the relative error of the produced interval (middle) and the spatial error of
the output. (l) Performance of the network during the training of t-DM,
quantified by the probability to successfully produce time interval (upper),
the relative error of the produced interval (middle) and the probability of
choice error (lower).  Figure S2: Interval production task. (a) Trajectory
speed with time in the perception epoch, shaded belt indicating s.e.m.
(standard error of mean). (b) Probability distribution function (p.d.f) of
scaling indexes of the activities of single neurons in the production epoch,
after counting neurons with the top 10% highest activity (upper panel), top
50% (middle panel) and all neurons (lower panel). (c) The scaling index and
explained variance of principal components (PC) in the production epoch. (d)
We calculated the scaling components in the subspace spanned by the first nine
principal components. Shown are the first (upper) and last (lower) scaling
component of the production epoch of an example training configuration. Color
of lines indicate to-be-produced interval $T$. (e) The mean activity of the
last scaling component as a function of $T$, with the activities when $T=600$
ms and $T=1200$ ms are respectively normalized to be 0 and 1. (f) Scaling
index (blue) and ratio of explained variance (orange) in the subspace spanned
by the accumulated scaling components. This panel is in the same style as Fig.
2n, except that it analyzes the perception epoch of IP task. (g,h) These two
panels explain the relationship between the low dimensionality of manifold
$\mathcal{M}$ at the end of the delay epoch and the dominance of neurons
monotonically tuned by $T$ (Section S2). (g) Histogram of the elements of PC1
of the manifold $\mathcal{M}$ at the end of the delay epoch at different $T$s
of an example training configuration. Note that the elements corresponding
with monotonically decreasing (MoD) and monotonically increasing (MoI) neurons
have different signs. (h) In 16 training configurations, for a given element
in PC1 of $\mathcal{M}$, it has over 98% probability to have the same sign
with most other elements corresponding with neurons of the same type, while
have the opposite sign with most other elements corresponding with neurons of
the opposite type. In panels c,e, error bars indicate s.e.m. over training
configurations. Figure S3: Interval comparison tasks. (a-c) Stimulus1 epoch.
(a) Population activity in the stimulus1 epoch in the subspace of the first
three PCs. Colors indicate the duration $T$ of the epoch. Stars and circles
respectively indicate the starting and ending points of the stimulus1 epoch.
(b) Coefficient of determination ($R^{2}$) that quantifies the overlap of the
firing profiles of individual neurons at different $T$s, in the same style as
Fig. 2d in the main text. (c) Trajectory speed as a function of time in the
stimulus1 epoch, shaded belt indicating s.e.m. (d-h) Delay epoch. (d)
Trajectory speed in the delay epoch when $T=600$ ms (blue) and 1200 ms (red),
in the same style as Fig. 2f. (e) Ratio of explained variance of the first
five PCs of manifold $\mathcal{M}$ at the end of the delay epoch, in the same
style as Fig. 2g. (f) The position of the state at the end of the delay epoch
projected in the first PC of manifold $\mathcal{M}$ as a function of $T$, in
the same style as Fig. 2h. (g) The distance between two adjacent curves in the
delay epoch as a function of time, in the same style as Fig. 2i. (h) The
portions of monotonically decreasing (MoD), monotonically increasing (MoI),
and non-monotonic (non-M) types of neurons at the end of the delay epoch, in
the same style as Fig. 2k. (j-o) Stimulus2 epoch. (i) Population activity in
the stimulus2 epoch in the subspace of the first three PCs. The meanings of
color scheme, stars and circles are the same as panel a. Triangles indicate
critical points. The duration of stimulus 2 is kept at 1200 ms. (j) Scaling
index (blue) and ratio of explained variance (orange) in the subspace spanned
by the accumulated scaling components, in the same style as Fig. 2n. In this
panel and panels k-n, only the trajectories from the beginning of stimulus 2
to the critical points are studied. (k) Trajectory speed in the subspace of
the first three scaling components, in the same style as Fig. 2o. (l)
Probability distribution of the scaling indexes of single neurons, in the same
style as Fig. S2b. (m) The scaling index and explained variance of principal
components, in the same style as Fig. S2c. (n) Mean activity of the last
scaling component as a function of $T$, in the same style as Fig. S2e. (o)
Left panel: speed of the trajectory before (blue) and after (red) the critical
point in the subspace of the first three scaling components (SC). SCs are
calculated using the trajectories before the cirtical points, the red line is
plotted by projecting the trajectories after the critical points into the
subspace of SCs calculated using those before critical points. Right panel:
speed of the trajectory before (blue) and after (red) the critical point in
the full population state space.  Figure S4: Timed spatial reproduction task.
(a,b) Perception epoch. (a) Coefficient of determination ($R^{2}$) that
quantifies the overlap of the firing profiles of individual neurons at
different $T$s in the perception epoch, in the same style as Fig. 2d. (b)
Trajectory speed as a function of time in the perception epoch, shaded belt
indicating s.e.m. (c-f) Delay epoch. (c) Trajectory speed as a function of
time in the delay epoch when $T=600$ ms (blue) and 1200 ms (red), in the same
style as Fig. 2f. (d) The manifold $\mathcal{M}$ at the end of the delay epoch
are parameterized by both time interval $T$ between the first two pulses and
the spatial location $x$ of the first pulse. We denote $\mathcal{M}(T;x_{0})$
(or $\mathcal{M}(x;T_{0})$) to be the set of dots in $\mathcal{M}$ at specific
location $x_{0}$ (or time interval $T_{0}$). Left panel: the position of the
state at the end of the delay epoch projected to the first PC of
$\mathcal{M}(T;x_{0})$ as a function of $T$, with the position when $T=600$ ms
(or 1200 ms) normalized to be 0 (or 1), in the same style as Fig. 2h. Gray
curves: results from 16 training configurations, each at a randomly chosen
$x_{0}$. Blue curve: mean value averaging over $x_{0}$ and training
configurations. Right panel: the position of the state in the first PC of
$\mathcal{M}(x;T_{0})$. We see that in most training configurations, the
position in $\mathcal{M}(x;T_{0})$ encodes $x$ continuously and linearly, but
big jump happens in some configurations. (e) The distance between two adjacent
curves in the delay epoch as a function of time, similar to Fig. 2i. (f) The
portions of monotonically decreasing (MoD), monotonically increasing (MoI) and
non-monotonic (non-M) types of neurons tuned by $T$ at the end of the delay
epoch, in the same style as Fig. 2k. (g-k) Production epoch. (g) Scaling index
(blue) and ratio of explained variance (orange) in the subspace spanned by the
accumulated scaling components in the production epoch, averaging over spatial
locations and training configurations, in the same style as Fig. 2n. (h)
Trajectory speed in the subspace of the first three scaling components in
production epoch, in the same style as Fig. 2o. (i) Probability distribution
of the scaling indexes of single neurons, in the same style as Fig. S2b. (j)
The scaling index and explained variance of principal components, similar to
Fig. S2c. (k) Mean activity of the last scaling component, similar to Fig.
S2e. Error bars representing s.e.m. are much smaller than the plot markers.
Figure S5: Timed decision making task. (a-c) Perception epoch. (a) Left:
Firing profiles of two example neurons in the perception epoch. Colors
indicate $c$ value, which is the half difference between the strength of the
presented stimuli. Right: Trajectories in the subspace of the first two PCs.
Stars and circles respectively indicate the starting and ending points of the
perception epoch. (b) Coefficient of determination ($R^{2}$) that quantifies
the overlap of the firing profiles of individual neurons at different $T$s in
the perception epoch, in the same style as Fig. 2d. (c) Trajectory speed as a
function of time in the perception epoch, shaded belt indicating s.e.m. (d-h)
Delay epoch. (d) Trajectories in the subspace of the first three PCs. Stars
and circles respectively indicate the starting and ending points of the delay
epoch. Blackness of circles indicates $T$ value as annotated. Curve color
indicates $c$ value as indicated in the color map of panel a, only $c=-0.04,$
-0.01, 0.01, 0.04 cases are plotted. (e) Trajectory speed as a function of
time in the delay epoch when $T=600$ ms (blue) and 1200 ms (red), in the same
style as Fig. 2f. (f) The position of the state in the first PC of
$\mathcal{M}(T;d_{0})$ as a function of $T$, with the position when $T=600$ ms
(or 1200 ms) normalized to be 0 (or 1), in the same style as Fig. 2h. Here,
$\mathcal{M}(T;d_{0})$ represents the set of dots in manifold $\mathcal{M}$ at
the end of the delay epoch at specific decision choice $d_{0}$. (g) The
distance between two adjacent curves in the delay epoch as a function of time,
in a similar style to Fig. 2i. Left panel: the two adjacent curves have the
same $c$ value, but slightly different $T$ values. Right panel: the two
adjacent curves have the same $T$ value, but different $c$ values. In the
right panel, blue (orange) curve represents the case when their $c$ values
have the same (different) sign, so that they have the same (different)
decision choice. We see that two trajectories representing the same
(different) choice tend to get close to (far away from) each other, consistent
with the scenario in panel d. (h) The portions of monotonically decreasing
(MoD), monotonically increasing (MoI) and non-monotonic (non-M) types of
neurons tuned by $T$ at the end of the delay epoch, in the same style as Fig.
2k. (i-m) Production epoch. (i) Scaling index (blue) and ratio of explained
variance (orange) in the subspace spanned by the accumulated scaling
components, averaging over $c$ values and training configurations, in the same
style as Fig. 2n. (j) Trajectory speed in the subspace of the first three
scaling components, in the same style as Fig. 2o. (k) Probability distribution
of the scaling indexes of single neurons, in the same style as Fig. S2b. (l)
The scaling index and explained variance of principal components, in the same
style as Fig. S2c. (m) Mean activity of the last scaling component, in the
same style as Fig. S2e. (n-s) The angle between first parameter-marginalized
principal components and mixed variances in the perception (panels n,o), delay
(panels p,q) and production epochs (panels r,s). These panels are in the same
style as Fig. 3d, e, g-j, except that the non-spatial information is decision
choice.  Figure S6: Examples that illustrate the geometry of the coding
combination of temporal and spatial information in t-SR. (a-e) Perception
epoch. (a) Each dot represents the angle between F-PC1 and S-PC1 as well as
their mixed variances in the perception epoch (after 400 ms of transient
period) of t-SR in a training configuration. (b-e) Iso-space lines in the
subspace spanned by F-PC1 and S-PC1, in the training configurations indicated
in panel a. Stars indicate the points after 400 ms of transient period from
the beginning of the perception epoch, and circles indicate the ending points
of the perception epoch. Redness from light to strong indicates the spatial
locations $x=0,2,4,\cdots,18$. (f-j) The same as panels a-e, except for
showing the iso-space lines in the manifold $\mathcal{M}$ at the end of the
delay epoch, in the subspace spanned by the first time-interval PC (I-PC1) and
S-PC1. Stars and circles indicate $T=600$ ms and 1200 ms cases respectively.
(k-o) The same as panels a-e, except that the iso-space lines in the
production epoch are shown. Stars indicate the points after 200 ms of
transient period from the beginning of the production epoch, and circles
indicate the ending points of the production epoch.  Figure S7: Decoding
generalizability in t-SR. (a-b) Perception epoch. (a) Upper: decoding error as
a function of $|x_{train}-x_{test}|$, after Decoder 1 (solid line) or Decoder
2 (dashed line) is trained to read the time elapsed from the beginning of the
perception epoch (i.e., temporal flow) using the state trajectory at spatial
location $x_{train}$, and then tested at spatial location $x_{test}$, in the
same style as Fig. 4g. Horizontal dashed line indicates chance level,
supposing the decorder works by random guess. Lower: The correlations between
the angle (AG) between the first temporal-flow PC and the first spatial PC,
the mixed variance (MV) between temporal flow and spatial information, the
error of Decoder 1 (DE1) and the error of Decoder 2 (DE2), in the same style
as Fig. 4d, f. Note that the correlation between AG and MV is approximately
zero, see Section S1.4.5 for this point. (b) Upper: Decoding error as a
function of $|t_{train}-t_{test}|$, after Decoder 1 (solid line) or Decoder 2
(dashed line) is trained to read the spatial location at time $t_{train}$
after the beginning of the perception epoch, and then tested at time
$t_{test}$. Lower: Correlations between AG, MV, DE1 and DE2. (c-d) Delay
epoch. (c) Similar to panel a, except for decoding time interval across
spatial information using the state in manifold $\mathcal{M}$ at the end of
the delay epoch. (d) Decoding spatial information across time interval using
the states in manifold $\mathcal{M}$ at the end of the delay epoch. (e-h)
Production epoch. (e) Decoding temporal flow across spatial information in the
production epoch. The decoder was trained using
$\mathbf{r}(t;x_{train},T_{0})$ and tested using
$\mathbf{r}(t;x_{test},T_{0})$, where $\mathbf{r}(t;x_{0},T_{0})$ represents
the population activity as a function of $t$ at specific spatial information
$x_{0}$ and time interval $T_{0}$. $T_{0}=1200$ ms in this panel and panels f.
(f) Decoding space across temporal flow in the production epoch. The decoder
was trained using $\mathbf{r}(x;t_{train},T_{0})$ and tested using
$\mathbf{r}(x;t_{test},T_{0})$, where $\mathbf{r}(x;t_{0},T_{0})$ represents
the population activity as a function of spatial information $x$ at specific
time point $t_{0}$ and time interval $T_{0}$. (g) Decoding temporal flow
across time interval in the production epoch. The decoder was trained using
$\mathbf{r}(t;T_{train},x_{0})$ and tested using
$\mathbf{r}(t;T_{test},x_{0})$. The results are averaged over
$x_{0}\in[0,20]$. Upper left: The decoded value $t_{dec}$ as a function of the
time $t$ elapsed from the beginning of the production epoch, after Decoder 1
(solid line) or Decoder 2 (dashed line) was trained to read $t$ at $T=1200$
ms, and then tested at $T=600$ ms (blue), 900 ms (red) and 1200ms (yellow).
The dashed line indicates perfect temporal scaling. Upper right: Decoding
error as a function of $T$, after a decoder is trained to read scaled temporal
flow $t/T$ at $T=1200$ ms (indicated by the vertical dashed line), and then
tested at $T=T_{1}$. Lower: correlations. (h) Decoding space across time
interval in the production epoch. The decoder was trained using
$\langle\mathbf{r}(x;T_{train},t_{0})\rangle_{t_{0}}$ and tested using
$\langle\mathbf{r}(x;T_{test},t_{0})\rangle_{t_{0}}$, where
$\langle\cdot\rangle_{t_{0}}$ means averaging over temporal flow $t_{0}$.
Figure S8: Decoding generalizability in t-DM. All panels are in the same style
as Fig. S7, except that the non-temporal information in t-DM is the decision
choice. Note that in some panels (lower panels of b, d, h), the correlation
between DE2 and AG as well as the correlation between DE2 and MV are absent.
The reason is that in these cases, the decoding error is perfectly zero in all
training configurations, so the correlation is undefined.  Figure S9:
Sequential activity and network structure. (a) The neuronal activity (with
maximum normalized to 1) in the production epoch of IP task in an example
training configuration, sorted according to peak time. (b, c) The same as
panel a, but for the stimulus1 (panel b) or stimulus2 (panel c) epoch of IC.
(d) Mean (solid line) and s.d. (shaded belt) of the recurrent weights as a
function of the peak order difference between post- and pre-synaptic neurons
in the production epoch of IP. (e, f) The same as panel d, but for the
stimulus1 (panel e) or stimulus2 (panel f) epoch of IC. (g) Recurrent weight
as a function of the difference $|x_{1}-x_{2}|$ between the preferred spatial
locations of post- and pre-synaptic neurons and their peak order difference in
the production epoch of t-SR. (h) Recurrent weight as a function of peak order
difference in the sequence of neurons with the same (blue) or different
(orange) preferred decision choices in the production epoch of t-DM. Shaded
belt indicates s.e.m.  Figure S10: Coding geometry and network structure in
the absence of timing task requirement. (a,b) The angle and mixed variance
between the subspaces coding temporal flow (F), time interval (I) and spatial
information (S) in the delay epoch of t-SR, in the same style as Fig. 3i, j.
(c,d) Similar to panel a,b, except for the delay epoch of t-DM, where the non-
temporal information is decision choice (D). (e) The angle between the first
temporal-flow PC and the first spatial (in t-SR, SR, COMP and CD) or decision-
choice (in t-DM, DM and cue-DM) PC. Whisker plots: center line, median; box,
25th to 75th percentiles; whiskers, $\pm 1.5\times$ the interquartile range.
In t-SR and t-DM, the perception epoch is studied; in SR, COMP and CD, the
delay epoch is studied; in DM and cue-DM, the stimulus-presentation epoch is
studied. Asterisk indicates significant ($p<0.05$) larger than $45^{\circ}$ (t
test). The horizontal dotted line indicates $45^{\circ}$, the vertical dotted
line separates the spatial task group (t-SR, SR, COMP and CD) from the
decision-making task group (t-DM, DM and cue-DM). The two horizontal dashed
line indicate the median values of t-SR and t-DM (which respectively are the
only timing task in each group) separately. (f) Mixed ratio $\rho$ in several
tasks, where $\rho=v_{\text{min}}/\min(v_{\text{time}},v_{\text{non-time}})$,
where $v_{\text{min}}$ is the mixed variance, $v_{\text{time}}$ and
$v_{\text{non-time}}$ are the variance explained by temporal and non-temporal
information separately. (g) Recurrent weight as a function of the difference
$|x_{1}-x_{2}|$ between the preferred spatial locations of post- and pre-
synaptic neurons and their peak order difference in the delay epoch of SR. (h)
The same as panel g, except for COMP. (i) The same as panel g, except for CD.
(j) Recurrent weight as a function of peak order difference in the sequence of
neurons with the same (blue) or different (orange) preferred decision choices
during the presentation of the stimuli in cue-DM. Shaded belt indicates s.e.m.
(k) The same as panel j, except for DM. Figure S11: Dynamics of the network
when trained to produce long time intervals. (a-b) Perception epoch. (a)
Population activity in the perception epoch in the subspace of the first three
PCs. Colors indicate the time interval $T$. Stars and circles respectively
indicate the starting and ending points of the perception epoch. (b)
Coefficient of determination ($R^{2}$) that quantifies the overlap of the
firing profiles of individual neurons at different $T$s in the perception
epoch, in the same style as Fig. 2d. (c-g) Delay epoch (c) Trajectory speed in
the delay epoch when $T=1200$ ms (blue) and 2400 ms (red), in the same style
as Fig. 2f. (d) Ratio of explained variance of the first five PCs of manifold
$\mathcal{M}$ at the end of the delay epoch, in the same style as Fig. 2g. (e)
The position of the state at the end of the delay epoch projected in the first
PC of manifold $\mathcal{M}$ as a function of $T$, in the same style as Fig.
2h. (f) The distance between two adjacent curves in the delay epoch as a
function of time, in the same style as Fig. 2i. (g) The portions of
monotonically decreasing (MoD), monotonically increasing (MoI) and non-
monotonic (non-M) types of neurons tuned by $T$ at the end of the delay epoch,
in the same style as Fig. 2k. (h-l) Production epoch. (h) Scaling index (blue)
and ratio of explained variance (orange) in the subspace spanned by the
accumulated scaling components, in the same style as Fig. 2n. (i) Trajectory
speed in the subspace of the first three scaling components, in the same style
as Fig. 2o. (j) Probability distribution of the scaling indexes of single
neurons, in the same style as Fig. S2b. (k) The scaling index and explained
variance of principal components, in the same style as Fig. S2c. (l) Mean
activity of the last scaling component as a function of $T$, in the same style
as Fig. S2e.

## References

  * [1] R. Hattori, K. V. Kuchibhotla, R. C. Froemke, and T. Komiyama, “Functions and dysfunctions of neocortical inhibitory neuron subtypes,” Nat. Neurosci., vol. 20, pp. 1199–1208, 2017.
  * [2] J. Wang, D. Narain, E. A. Hosseini, and M. Jazayeri, “Flexible timing by temporal scaling of cortical responses,” Nat. Neurosci., vol. 21, pp. 102–110, 2018.
  * [3] A. E. Orhan and W. J. Ma, “A diverse range of factors affect the nature of neural representations underlying short-term memory,” Nat. Neurosci., vol. 22, pp. 275–283, 2019.
  * [4] G. R. Yang, M. R. Joglekar, H. F. Song, W. T. Newsome, and X.-J. Wang, “Task representations in neural networks trained to perform many cognitive tasks,” Nat. Neurosci., vol. 22, pp. 297–306, 2019.
  * [5] Q. V. Le, N. Jaitly, and G. E. Hinton, “A simple way to initialize recurrent networks of rectifed linear units,” arXiv https://arxiv.org/abs/1504.00941, 2015.
  * [6] A. Goel and D. V. Buonomano, “Temporal interval learning in cortical cultures is encoded in intrinsic network dynamics,” Neuron, vol. 91, pp. 1–8, 2016\.
  * [7] C. J. MacDonald, K. Q. Lepage, U. T. Eden, and H. Eichenbaum, “Hippocampal ’time cells’ bridge the gap in memory for discontiguous events,” Neuron, vol. 71, pp. 737–749, 2011.
  * [8] R. H. R. Hahnloser, A. A. Kozhevnikov, and M. S. Fee, “An ultra-sparse code underlies the generation of neural sequence in a songbird,” Nature, vol. 419, pp. 65–70, 2002.
  * [9] E. Pastalkova, V. Itskov, A. Amarasingham, and G. Buzsáki, “Internally generated cell assembly sequences in the rat hippocampus,” Science, vol. 321, pp. 1322–1327, 2008.
  * [10] C. D. Harvey, P. Coen, and D. W. Tank, “Choice-specifc sequences in parietal cortex during a virtual-navigation decision task,” Nature, vol. 484, pp. 62–68, 2012.
  * [11] D. Kobak, W. Brendel, C. Constantinidis, C. E. Feierstein, A. Kepecs, Z. F. Mainen, X.-L. Qi, R. Romo, N. Uchida, and C. K. Machens, “Demixed principal component analysis of neural population data,” eLife, vol. 5, p. e10989, 2016.
  * [12] J. D. Murray, A. Bernacchia, N. A. Roy, C. Constantinidis, R. Romo, and X.-J. Wang, “Stable population coding for working memory coexists with heterogeneous neural dynamics in prefrontal cortex,” Proc. Natl. Acad. Sci. USA, vol. 114, pp. 394–399, 2017.


Send correspondence to<EMAIL_ADDRESS>
# Instructive artificial intelligence (AI) for human training, assistance, and
explainability

Nicholas Kantack  The Johns Hopkins Applied Physics Laboratory, Laurel, MD
20723  Nina Cohen  The Johns Hopkins Applied Physics Laboratory, Laurel, MD
20723  Nathan Bos  The Johns Hopkins Applied Physics Laboratory, Laurel, MD
20723  Corey Lowman  The Johns Hopkins Applied Physics Laboratory, Laurel, MD
20723  James Everett  The Johns Hopkins Applied Physics Laboratory, Laurel, MD
20723  Timothy Endres  The Johns Hopkins Applied Physics Laboratory, Laurel,
MD 20723

###### Abstract

We propose a novel approach to explainable AI (XAI) based on the concept of
“instruction” from neural networks. In this case study, we demonstrate how a
superhuman neural network might instruct human trainees as an alternative to
traditional approaches to XAI. Specifically, an AI examines human actions and
calculates variations on the human strategy that lead to better performance.
Experiments with a JHU/APL-developed AI player for the cooperative card game
Hanabi suggest this technique makes unique contributions to explainability
while improving human performance. One area of focus for Instructive AI is in
the significant discrepancies that can arise between a human’s actual strategy
and the strategy they profess to use. This inaccurate self-assessment presents
a barrier for XAI, since explanations of an AI’s strategy may not be properly
understood or implemented by human recipients. We have developed and are
testing a novel, Instructive AI approach that estimates human strategy by
observing human actions. With neural networks, this allows a direct
calculation of the changes in weights needed to improve the human strategy to
better emulate a more successful AI. Subjected to constraints (e.g. sparsity)
these weight changes can be interpreted as recommended changes to human
strategy (e.g. “value A more, and value B less”). Instruction from AI such as
this functions both to help humans perform better at tasks, but also to better
understand, anticipate, and correct the actions of an AI. Results will be
presented on AI instruction’s ability to improve human decision-making and
human-AI teaming in Hanabi.

###### keywords:

XAI, explainability, interpretability, hanabi, human-machine, teaming,
instructive, instruction

##  1 INTRODUCTION

AI systems have demonstrated the ability to perform tasks remarkably well from
games [1, 2] to medical diagnosis [3]. Many of these AI are comprised of deep
neural networks from which it is very hard to extract insight or explanations
of the decision the network makes. Therefore, in many cases a neural network
can discover novel insights into a domain but cannot communicate these
insights to the humans that developed the network. This fundamental problem
has sparked the active field of research into explainable AI (XAI)[4], AI for
which there are some measures in place to facilitate human understanding of
the AI’s decision.

In some cases, the unexplainability of AI is a barrier to its use. Such cases
are those in which humans are agents who must collaborate with the AI (which
typically requires some level of common understanding) and those in which
humans are significant stakeholders (e.g. when the AI is recommending medical
treatment). A research effort at JHU/APL entitled “Learning to Read Minds”
studied this challenge in the context of human-machine teaming in the
collaborative card game Hanabi [5]. Hanabi is a sort of cooperative solitaire
with imperfect information that requires players (human or machine) to be able
to infer the knowledge, intentions, and future actions from the behavior of
their teammates [6]. Hanabi is a game for which the traditional process of
self-play optimization (i.e. training an AI through millions of games played
between copies of the same AI) does not lead to successful human-machine
performance [5], primarily because AI agents can develop obscure conventions
(e.g. repurposing an in-game clue to mean something entirely different from
its semantic meaning) that will be automatically understood by their mirror
image during self-play, but completely incomprehensible to a human. This is
why agents such as the Simplified Action Decoder [7], Rainbow [8], and
Fireflower [9] often achieve perfect scores in self-play, yet achieve low
scores when playing with human teammates [5]. Furthermore, due to the lack of
effective XAI techniques, there appear no practical means for these complex
self-play AIs to explain these obscure conventions to humans (setting aside
whether humans are even capable of implementing these conventions once
understood).

The “Learning to Read Minds” research project included a JHU/APL-internal
challenge tasking staff with developing AI agents that would excel when
playing Hanabi with human strangers. The winning JHU/APL agent not only
achieved human-play scores higher than any found in literature to date, [10,
11, 12] but it did so in a way that was constrained to allow human-readable
descriptions of strategy (Figure 1). In particular, the JHU/APL agent
demonstrated the ability to develop deep insights into human strategy through
observation of human play, to understand how the human strategy interacted
with the agent’s strategy, and to adapt to discover a play style which
complements the human strategy. This study summarizes the agent’s structure
which enabled it to successfully collaborate with human teammates, and
introduces a novel type of explanation (we call “instruction”) to share AI
insights with human observers.

##  2 A Human-like Hanabi Agent

The JHU/APL agent (henceforth referred to as “agent”) was developed under the
philosophy that if it could play like humans, it would play well with humans.
The agent was designed to convert the input space of the game state to a
latent space of a small number of human-preferred factors (HPFs) which are
aspects of the game that humans are known to attend to when making decisions.
The agent utilizes twelve HPFs (Table 1) which were suggested by intermediate
Hanabi players. Constraining the attention of an AI in this fashion in order
to guarantee some level of interpretability after training is a known practice
[13, 14]. In the case of the JHU/APL Hanabi agent, an expected reward for each
possible action is computed based on the expected effect the action will have
on the HPFs. In particular, the expected value of an action is the inner
product of a factor vector $\vec{h}$ with a weights vector $\vec{w}$.
Therefore,

$\displaystyle y_{i}=\vec{h}_{i}^{t}\vec{w}\hskip 28.45274pt\implies\hskip
28.45274pt\vec{y}=H^{T}\vec{w}$ (1)

where $y_{i}$ is the expected reward for action $i$, and vectors $\vec{h}_{i}$
form the columns of $H$. The elements of $\vec{h}$ are the expected changes
that an action will induce on each of the HPFs (e.g. for the HPF of playing a
playable card, the corresponding element in $\vec{h}_{i}$ is the probability
that action $i$ will result in the playing of a playable card). The elements
of $\vec{w}$ are the relative values of each HPF with respect to one another.
Thus, while $H$ represents information about the game state, $\vec{w}$
represents the agent’s strategy. Altering the elements of $\vec{w}$ can
dramatically alter the play style of the agent.

On each move, the agent calculates $\vec{y}$ which stores the expected reward
for each possible action. The agent always chooses the action with the highest
expected reward among the legal actions available. Of note, this technique
does not involve and consideration of moves beyond the ply under
consideration. Rather, the agent is pursuing an immediate improvement of the
game state with respect to the chosen HPFs.

Table 1: The twelve factors, along with their values for three different JHU/APL agent play styles (human-like, human-complementary, and self-play). To help interpret some of the factors, consider the following definitions. endangered card \- a card for which no copy has been played, yet there is only one copy of this card remaining in play. unneeded card \- A card which cannot be played in the future, for any reason. Factor | Weights
---|---
human-like | human-compl. | self-play
Playing a playable card | 1 | $\infty$ | 11
Playing unplayable card |  |  |
(fewer than 2 strikes) | -1 | -1 | -1
Playing an unplayable card (2 strikes) | 3 | $\infty$ | 3
Other player playing a |  |  |
playable card | 1.5 | 10 | 2
Other player playing an |  |  |
unplayable card | 0 | 0 | 1
Discarding a non-endangered card | 0.1 | 0.55 | 0.8
Discarding an unneeded card | 0.25 | 1 | 0
Playing a singled out card | 3 | 1.5 | 5
Giving a clue that singles |  |  |
out a playable card | 3 | 3 | 2
Giving a clue that singles |  |  |
out a non-playable card | 0 | -5 | 4
Discarding a singled out card | -0.5 | -2 | -3
Added value to any clue |  |  |
per info token held | 0.5 | 0.1 | 0
Figure 1: The score distributions are shown for Simplified Action Decoder (a
special off-belief version made for the competition), Rainbow, Fireflower, and
the JHU/APL agent. These scores were obtained by pairing an agent with a human
teammate (drawn from a pool of 21).

###  2.1 Modeling Human Decision Making

While human-like play was the preliminary goal during the agent’s development,
the first training efforts were aimed at generating decent self-play scores.
For this phase, the training of the agent was separated into epochs. Each
epoch consists of a four-dimensional, full factorial design experiment on a
subset of four elements from $\vec{w}$. Each element under test was given
three test values (a low, medium, and high) around the neighborhood of where
the optimal value was expected to be. Therefore, each epoch tested $3^{4}=81$
unique $\vec{w}$ vectors. For each $\vec{w}$ vector tested, 200 games were
played between identical copies of the agent. The elements of $\vec{w}$ under
test were not altered until a an epoch occurred for which the highest score
was achieved by assigning the medium value for each element under test (i.e.
increasing or decreasing any element led to poorer performance). The
progression of self-play scores during this development phase are shown in
Figure 2.

Once the agent was optimized for self-play, the next objective was to find a
strategy vector $\vec{w}$ that would lead to play that was as human-like as
possible. To facilitate this exploration, a dataset of 376 decisions was
collected by examining the play of one of the authors. With a dataset of
decision made by a single human, we aimed to determine if a strategy vector
$\vec{w}$ could be fitted to a particular human’s play style rather than a
$\vec{w}$ which represented some ambiguous (perhaps bad) play style that was
averaged across humans with potentially dissimilar play styles. The
progression of increasing humanness is displayed in Figure 3. The highest
humanness fraction of any agent was 64.2%, achieved by the human-like agent
(that is, the agent was able to independently agree with the human decision in
64.2% of the game states examined).

Once the human-like version of the agent had been optimized for fitting the
dataset of human decisions, a final training effort was made by pairing a
training version of the agent with the human-like version. In this fashion,
the training process was intended to approximate playing with a human
teammate. As before, full factorial design experiments were run altering four
elements of $\vec{w}$ per epoch, each across three levels. At then end of the
training process, the “human-complementary” version of the agent was created.
Cross play results (Figure 4) illustrate the performance of different
combinations of agents developed.

###  2.2 An Important Note on Human Perception of Strategy

It is worth noting that the JHU/APL agent needed to make significant changes
to its initial $\vec{w}$ vector in order to accurately predict human decision
making in the game, despite the fact that the initial $\vec{w}$ given to the
agent was intended to accurately describe human decision making. For this
reason, it became clear that human players could not accurately depict the
weights they attributed to HPFs. This has profound implications for XAI. A
common XAI approach would have involved taking the values from the self-play
strategy vector in Table 1 and describing these to a human player (e.g. “You
should value discarding a non-endangered card at 0.8”). However, if a human
already egregiously misunderstands what value they actually attribute to these
HPFs, it is unlikely that the human will be able to act on this insight.
Rather, it would perhaps be more suitable for us to look at the difference in
weights between the human-like agent and the self-play agent, since doing so
would allow us to specify corrections a human should make to their strategy
(e.g. “you should value discarding a non-endangered card more”). These
corrections are human interpretable regardless of whether the human accurately
understands their current strategy. This is the principal idea behind AI
instruction.

Figure 2: Self play scores are shown for different versions of the agent
during the initial development process. Each transition from one version to
the next was accompanied by changes to the weights vector $\vec{w}$ or the
addition of new elements to the weights vector. Notably, the human-like agent
had significantly poorer self-play scores, consistent with the fact that the
human-like agent was optimized to agree with a database of human decisions.
Figure 3: The humanness of major versions of the agent are shown. Humanness is
defined as the fraction of human decisions with which the agent agrees when
analyzing a database of 376 game decisions made by one of the authors. Of the
agents shown, only the “human-like” agent was explicitly optimized for maximal
humanness. The considerable humanness of the other models indicates how
optimizing an HPF focused agent for self play can lead to considerably human-
like performance. Figure 4: Average scores are shown for different pairings of
the agents (and humans). The human-human score is included for reference, but
is an average across a small number of games played within the development
team. All other scores are averages across at least 10 games. Error margins
for each bin are less than 1 point. It is clear that the human-complementary
agent achieves the highest score of any agent paired with the human-like agent
(since these are the exact conditions under which the human-complementary
agent was optimized). It is notable that the average human + self-play and
human + human-compl. scores are identical in this plot (where a difference is
expected), but it is also worth noting that a very small number (two) of very
experienced humans were represented in these data. Therefore, while these
scores are useful for comparing how non-human agents perform in different
pairings, generalizations about human play from this figure should be made
with great caution.

##  3 Theory of AI Instruction

AI instruction is defined in the context of explaining differences in strategy
in the form of changes on weights. Therefore, it is relevant to consider a
difference in outputs (say, $\vec{y}_{1}$ from $\vec{w}_{1}$ and $\vec{y}_{2}$
from $\vec{w}_{2}$).

$\displaystyle\delta\vec{y}=\vec{y}_{1}-\vec{y}_{2}=H^{T}\vec{w}_{1}-H^{T}\vec{w}_{2}$
(2)
$\displaystyle\delta\vec{y}=H^{T}\left(\vec{w}_{1}-\vec{w}_{2}\right)=H^{T}\delta\vec{w}$
(3)

###  3.1 A note on Strategy vs. Perception

It is possible to imagine a difference in performance, $\delta\vec{y}$,
arising not from a difference in strategy ($\delta\vec{w}$), but rather a
difference in perception of the game state $\delta H$. This is particularly
practical if the elements of $H$ concern complex changes in the game state
such as the probabilities of certain outcomes (as it does in Hanabi). In this
case,

$\displaystyle\delta\vec{y}=\delta H^{T}\vec{w}$ (4)

In fact, there is ambiguity between this and (3), since $HH^{T}$ is a
$12\times 12$ matrix that will tend to be full rank, and thus the relation

$\displaystyle\delta\vec{w}=\left(HH^{T}\right)^{-1}H\delta H^{T}\vec{w}$ (5)

indicates that any strategic difference $\delta\vec{w}$ could be interpreted
instead as an observation error $\delta H$. This illustrates yet another
reason for providing instruction in the form of $\delta\vec{w}$ rather than an
explanation in the form of $\vec{w}_{2}$. By the very nature of this
equivalence relation, a recommended change in strategy, $\delta\vec{w}$ can
compensate for both misperception and strategic deficiency (and mixtures
thereof).

###  3.2 Non-uniqueness of and Constraints on $\delta\vec{w}$

Since $H^{T}$ is a tall matrix ($20\times 12$), $H^{T}$ has a non-empty null
space. Therefore, the condition

$\displaystyle\delta\vec{y}=H^{T}\left(\delta\vec{w}+\vec{n}\right)$ (6)

is satisfied for any $\vec{n}$ in the null space of $H^{T}$, and so any
$\delta\vec{w}+\vec{n}$ is a valid description of a strategic difference
needed to elicit the decision difference $\delta\vec{y}$. This non-uniqueness
of $\delta\vec{w}$ is advantageous, because it allows multiple possible
$\delta\vec{w}$ to be compared for fitness according to human friendly
constraints (e.g. norm minimality, sparsity, etc.).

###  3.3 Generating AI Instruction

Suppose that a human subject is presented with $g$ game states
$H_{1},...,H_{g}$, and that these game state matrices are stacked into a
$12\times 20\times g$ tensor $\mathcal{H}$. As before, a strategy $\vec{w}$
can be combined with a game state to yield a vector of outputs, $\vec{y}$.

$\displaystyle\vec{y}_{k}=\mathcal{H}(:,:,k)^{T}\vec{w}$ (7)

In the following formalism, a subscript $h$ corresponds to a human, while a
subscript $i$ corresponds to an ideal (typically a successful AI). Let’s
assume that the index of the maximum element of $\vec{y}_{h}$ indicates the
decision that will be taken (per strategy $\vec{w}_{h}$) for the game state
used. In this case, two different $\vec{y}$ vectors may still specify the same
action if their maximal elements occupy the same index. If not, then it is
worth describing the nearest (min $|\vec{y}_{i}-\vec{z}|$) vector $\vec{z}$
such that $\vec{y}_{h}$ and $\vec{z}$ have maximal elements in the same index
position (a position different than the maximal element of $\vec{y}_{h}$).
Suppose $y_{h}(t)=$max$(\vec{y}_{h})$, but no other information is known about
$\vec{y}_{h}$. Then let $m$ be the average of all the terms of $\vec{y}_{i}$
that are greater than $y_{i}(t)$. Then $\vec{z}$ is defined as

$\displaystyle z(i)=\left\\{\begin{matrix}y_{i}(j)&y_{i}(j)<m\ \&\ i\neq t\\\
m+\varepsilon&y_{i}(j)<m\ \&\ j=t\\\ m&y_{i}(j)>m\\\ \end{matrix}\right.$ (8)

where $\varepsilon$ is some small, positive tie-breaking factor. If the
vectors $\vec{z}$ (of which there are $g$) are made the columns of a $20\times
g$ matrix $Z$, and the output vectors $y_{i}$ from strategy $\vec{w}_{i}$ are
made columns of a $20\times g$ matrix $Y_{i}$, then we can relate $Z$ and the
game state tensor $\mathcal{H}$ as follows.

$\displaystyle\mathcal{H}(:,:,k)^{T}\delta\vec{w}=Z(:,k)-Y_{i}(:,k)\ \forall\
k\in[1,g]$ (9)

Where $\delta\vec{w}$ is a strategy change needed so that
$\vec{w}_{i}+\delta\vec{w}$ and $\vec{w}_{h}$ arrive at the same decision for
every game state in $\mathcal{H}$. To calculate for $\delta\vec{w}$, we can
utilize the following matrix unfolding.

$\displaystyle\left[\begin{matrix}\mathcal{H}(:,:,1)^{T}\\\
\mathcal{H}(:,:,2)^{T}\\\ ...\\\
\mathcal{H}(:,:,g)^{T}\end{matrix}\right]\delta\vec{w}=\widetilde{H}\delta\vec{w}=\text{vec}\left(Z-Y_{i}\right)$
(10)

This is an overspecified linear system, so a least squared error solution can
be taken for $\delta\vec{w}$. Then, $\delta\vec{w}$ is the norm-minimal change
to apply to $\vec{w}_{i}$ to better concur with strategy $\vec{w}_{h}$ in each
of the $g$ game states. If $\vec{w}_{h}$ is the strategy of a human, and
$\vec{w}_{i}$ is an ideal, $\delta\vec{w}$ is the change to the ideal needed
to concur with the human. The inverse ($-\delta\vec{w}$) has elements which
comprise the instructions that should be given to the human. In essence, the
instructed changes are the opposite of those needed for the ideal to be
altered to make the same decisions the human made.

###  3.4 Properties of the Generated $\delta\vec{w}$

$\delta\vec{w}$ is not guaranteed to produce consensus between the starting
strategy and the ideal when adopted. Formally, it does not always hold that

$\displaystyle\mathcal{H}(:,:,k)(\vec{w}_{i}+\delta\vec{w})=\mathcal{H}(:,:,k)\vec{w}_{h}$
(11)

However, it is possible (and desired) for this relation to circumstantially
hold for many $k$ values. Increasing $\varepsilon$ in the above formulation
will tend to increase the number of game states in which consensus is built
but at the expense of a larger norm $\delta{w}$ (i.e. bigger recommended
changes to $\vec{w}_{h}$). Even so, total consensus between the ideal and
modified strategies is rarely achieved because the model for decisions
generated from game states given by (1) may not accurately describe all
decisions ($\vec{y}$) made by a human (e.g. due to momentary misperception,
distraction, and attention to factors not captured in $H$). The quantity
$\lambda$ defined as

$\displaystyle\lambda=\min_{\vec{w}}\frac{1}{g}\sum_{k=1}^{g}|\vec{y}_{h}-\mathcal{H}(:,:,k)^{T}\vec{w}|$
(12)

may be introduced as a figure of merit for the list of factors which define
the strategy vector $\vec{w}$. Furthermore, $\lambda$ can be used to measure
the utility of elements of $\vec{w}$ by examining the change in $\lambda$
induced by the removal or inclusion of factors. Ideally, the only factors kept
would be those whose inclusion result in a significant decrease in $\lambda$.

Similarly, one can define a figure of merit for generated instruction. If we
define $f(\vec{a},\vec{b})$ as

$\displaystyle
f(\vec{a},j)=\left\\{\begin{matrix}1&a(j)=\text{max}(\vec{a})\\\
0&\text{otherwise}\end{matrix}\right.$ (13)

If $n_{k}$ is the index of the decision the human instructee made for game
state $k$, then it is possible to evaluate the quality $q(\delta\vec{w})$ of
instructions as

$\displaystyle
q(\delta\vec{w})=\frac{1}{g}\sum_{k=1}^{g}f(\mathcal{H}(:,:,k)^{T}(\vec{w}_{i}+\delta\vec{w}),n_{k})$
(14)

$q(\delta\vec{w})$ falls in $[0,1]$ and can be interpreted as the fraction of
human decisions that can be understood as a variation ($\delta\vec{w}$) on an
ideal ($\vec{w}_{i}$).

###  3.5 Full AI Instruction Algorithm with Quality Monitoring

We recommend the algorithm in Figure 5 for generating AI instruction. The
algorithm has two preparation steps. The first is to train up an ideal
strategy ($\vec{w}_{i}$), and the second is to aggregate a dataset of human
decisions $\vec{n}$ paired with the game states in which they were made (slabs
of the $\mathcal{H}$ tensor). While it is possible to terminate the algorithm
after the step that assigns $\delta\vec{w}$, this algorithm includes a post-
processing component which seeks to zero out as many elements of
$\delta\vec{w}$ as possible while maintaining some preset explanatory fidelity
$\alpha$ to the human decision set. The purpose of this post-processing is to
generate instruction which concerns changes in as few of values as possible.
This is motivated by the assumption that low dimensional instructions are
easier for humans to understand (i.e. require focusing on fewer aspects of the
game in subsequent play).

Instructive AI Algorithm

$\alpha\leftarrow\text{Chose accuracy threshold in [0,1]}$

$\vec{w}_{i}\leftarrow\text{Train}(\hat{\mathcal{H}},C)$ Learn ideal weights
on a dataset $\hat{\mathcal{H}}$

$\mathcal{H},\vec{n}\leftarrow\text{Observe }g\text{ human decisions}$

for $k\in[1,g]$ do

$Y_{i}(:,k)\leftarrow\mathcal{H}(:,:,k)^{T}\vec{w}_{i}$

$Z(:,k)\leftarrow(Y_{i}(k),n(k))$ per (8)

end for

$\delta\vec{w}\leftarrow\text{Solve
}\widetilde{H}\delta\vec{w}=\text{vec}\left(Y_{i}-Z\right)$

$q\leftarrow\frac{1}{g}\sum_{k=1}^{g}f(\mathcal{H}(:,:,k)^{T}(\delta\vec{w}_{i}+\delta\vec{w}),n_{k})$

while $q>\alpha$ do

$\delta\vec{w}\leftarrow$ zero out element with smallest impact on $q$

$q\leftarrow\frac{1}{g}\sum_{k=1}^{g}f(\mathcal{H}(:,:,k)^{T}(\delta\vec{w}_{i}+\delta\vec{w}),n_{k})$

end while

Give $-\delta\vec{w}$ to human

Figure 5: This algorithm generates sparse AI instruction tailored to a set of
human decisions.

##  4 Experimental Results

During the “Learning to Read Minds” challenge, a database of 376 human
decisions in Hanabi games was generated. Preliminary results are shown based
on analysis of this dataset. To illustrate the instruction generation process,
a trial agent was created by copying the self-play agent. Because the self-
play agent already agrees with human decision at a high rate, the weight for
the non-endangered discard was inflated (to a value of 10). Then, in an
iterative process, instructions were generated (on how the trial agent could
better emulate human decision making based on the dataset), the trial agent
applied the instructed changes to its weights, and a new set of instructions
were generated. This process is shown to lead to asymptotic improvement in the
agreement between the trial agent and the human dataset (Figure 6).

High agreement ($68\%$) was achievable after 40 instruction based weight
updates. Furthermore, the spurious discard weight was shown to be brought into
closer agreement with the target strategy. Importantly, this (and other
initially matching weights) were shown to drift to new equilibrium values.
This serves as an empirical demonstration of the non-uniqueness of strategies
as described in the previous section. However, it is important to note that
the generation of a norm minimally different $Z$ matrix may not provide a
linear system in (10) that admits a solution that produces high prediction
accuracy when observing the target strategy. This is because the matrix $Z$
may be a poor estimation for the target strategy’s output vectors, a
circumstance that is increasingly likely when the instructee strategy differs
significantly from the target strategy.

These results (Figure 6) indicate that AI instruction can indeed provide
stepwise improvements to strategy which, taken iteratively, can lead to
significant improvement in the agreement between the instructee strategy and
the ideal. In this way, instructions serve as something of a proxy gradient of
a cost function, namely, agreement with the ideal. Utilizing the instructions
as a gradient for agent training was shown in this experiment to lead to
better humanness scores ($68\%$ vs. $64\%$) in a much shorter computation time
(minutes vs. hours) compared to the full factorial approach. Additionally,
these instructions provide a novel approach to portraying AI insight in a way
that is understandable to human observers. Specifically, these instructions
can be phrased as corrections to the weights attributed to human-preferred
factors, allowing for AI systems to develop an understanding of human decision
making and to share those insights through tailored instructions.

Figure 6: The effect of following successive batches of AI instructions are
shown. For this experiment, a trial version of the self-play agent was created
that had a non-endangered discard weight of 10 rather than the normal 0.8.
Instructions are generated for how this agent can better emulate the self-play
agent. (a) The agreement between the agents is shown as a function of how many
batches of instructions were generated. After each batch, the trial agent
applies the correction to its weights and reexamines agreement. (b) The trial
agent’s inflated discard weight decreases with time. Notably, it does not
appear to stabilize at the same value as the ideal (0.8). (c) The weight for
discarding unneeded cards is shown over the experiment. Note that while the
two agents initially had the same value (0), this weight drifts to a new value
as a result of the instructions, indicating that a new set of weights is being
discovered that agrees with the self-play agent over the decisions studied.
(d) Another weight is shown (which does not pertain to discards) to further
illustrate how instruction may not encourage the same, unique weights as the
ideal.

##  5 Conclusion

Leveraging insights obtained from the development of a highly successful,
artificially intelligent human teammate for Hanabi, we propose a technique of
instructive AI to better enable humans to obtain insight from complicated AI
systems. There are assumptions in this approach that may not hold true for
certain contexts. For instance, this technique hopes that the requisite
$\delta\vec{w}$ for consensus building is small. If not, then implementing a
$\delta\vec{w}$ may be just as confusing for humans as being told
$\vec{w}_{i}$, or perhaps even more so. More fundamental, the model given by
(1) may not accurately capture a majority of a human’s decisions, and
$\vec{w}$ is always at risk of missing elements that are crucial to a human’s
decision making. In general, it is challenging to produce a complete set of
values relevant to human decision making. For the purposes of this experiment,
the list of values is produced from human introspection and trial and error.
Techniques to organically learn the needed values may be possible and highly
valuable to the task of generating AI instruction, but are beyond the scope of
the experiments described above.

Many of the challenges described above apply in similar form to other methods
of XAI. However, instructive AI shows promise to circumvent some of the
greatest challenges of XAI and provide a novel framework in which further
research might push the frontier on extracting human-useful insight from
complex AI systems.

## References

  * [1] Silver, D., Schrittwieser, J., and et al., K. S., “Mastering the game of go without human knowledge,” Nature 550, 354–359 (2017).
  * [2] Silver, D., Hubert, T., Schrittwieser, J., Antonoglou, I., Lai, M., Guez, A., Lanctot, M., Sifre, L., Kumaran, D., Graepel, T., Lillicrap, T., Simonyan, K., and Hassabis, D., “Mastering chess and shogi by self-play with a general reinforcement learning algorithm,” (2017).
  * [3] McKinney, S. M. and et al., M. S., “International evaluation of an ai system for breast cancer screening,” Nature 577, 89–94 (2020).
  * [4] Guidotti, R., Monreale, A., Ruggieri, S., Turini, F., Pedreschi, D., and Giannotti, F., “A survey of methods for explaining black box models,” (2018).
  * [5] Anonymous, “Towards human-compatible ai teammates,” (2022).
  * [6] Bard, N., Foerster, J. N., Chandar, S., Burch, N., Lanctot, M., Song, H. F., Parisotto, E., Dumoulin, V., Moitra, S., Hughes, E., and et al., “The hanabi challenge: A new frontier for ai research,” Artificial Intelligence 280, 103216 (Mar 2020).
  * [7] Hu, H. and Foerster, J. N., “Simplified action decoder for deep multi-agent reinforcement learning,” (2021).
  * [8] Hessel, M., Modayil, J., van Hasselt, H., Schaul, T., Ostrovski, G., Dabney, W., Horgan, D., Piot, B., Azar, M., and Silver, D., “Rainbow: Combining improvements in deep reinforcement learning,” (2017).
  * [9] “Fireflower.” https://github.com/lightvector/fireflower (2018).
  * [10] Hu, H., Lerer, A., Peysakhovich, A., and Foerster, J., “”other-play” for zero-shot coordination,” (2021).
  * [11] Eger, M., Martens, C., and Cordoba, M. A., “An intentional ai for hanabi,” in [2017 IEEE Conference on Computational Intelligence and Games (CIG) ], 68–75 (2017).
  * [12] Siu, H. C., Pena, J. D., Chen, E., Zhou, Y., Lopez, V. J., Palko, K., Chang, K. C., and Allen, R. E., “Evaluation of human-ai teams for learned and rule-based agents in hanabi,” (2021).
  * [13] Yang, Z., Zhang, A., and Sudjianto, A., “Enhancing explainability of neural networks through architecture constraints,” (2019).
  * [14] Alvarez-Melis, D. and Jaakkola, T. S., “Towards robust interpretability with self-explaining neural networks,” (2018).