††thanks: These authors contributed equally††thanks: These authors contributed
equally††thanks: These authors contributed equally

# Analytical prediction of specific spatiotemporal patterns in nonlinear
oscillator networks with distance-dependent time delays

Roberto C. Budzinski  Tung T. Nguyen  Gabriel B. Benigno  Department of
Mathematics, Western University, London, ON, Canada  Brain and Mind Institute,
Western University, London, ON, Canada  Western Academy for Advanced Research,
Western University, London, ON, Canada  Jacqueline Đoàn  Department of
Mathematics, Western University, London, ON, Canada  Brain and Mind Institute,
Western University, London, ON, Canada  Western Academy for Advanced Research,
Western University, London, ON, Canada  Ján Mináč  Department of Mathematics,
Western University, London, ON, Canada  Western Academy for Advanced Research,
Western University, London, ON, Canada  Terrence J. Sejnowski  The Salk
Institute for Biological Studies, La Jolla, CA, USA  Division of Biological
Sciences, University of California at San Diego, La Jolla, CA, USA  Lyle E.
Muller<EMAIL_ADDRESS>Department of Mathematics, Western University, London,
ON, Canada  Brain and Mind Institute, Western University, London, ON, Canada
Western Academy for Advanced Research, Western University, London, ON, Canada

###### Abstract

We introduce an analytical approach that allows predictions and mechanistic
insights into the dynamics of nonlinear oscillator networks with heterogeneous
time delays. We demonstrate that time delays shape the spectrum of a matrix
associated to the system, leading to the emergence of waves with a preferred
direction. We then create analytical predictions for the specific
spatiotemporal patterns observed in individual simulations of time-delayed
Kuramoto networks. This approach generalizes to systems with heterogeneous
time delays at finite scales, which permits the study of spatiotemporal
dynamics in a broad range of applications.

##  I Introduction

What is the effect of heterogeneous time delays in networked systems? This
question is difficult to treat analytically in the context of multiple
distributed time delays. In recent work [1], we studied intracranial
electrophysiological recordings from human clinical patients during sleep. We
found that the 11-15 Hz sleep “spindle” oscillation, a brain rhythm important
for learning and memory [2], was not perfectly synchronized with zero phase
difference across the cortex; rather, sleep spindles are organized into
rotating waves that travel in a preferred direction (see Movie 1 in [1]).
Importantly, the propagation speed of the observed waves is consistent with
the axonal conduction speed of the long-range fiber network in the cortex (3-5
ms-1 [3]). This set of observations raises an important question: how do these
fibers, with no major anisotropy, create a specific spatiotemporal structure
with a preferred chirality?

In this work, we analyze a time-delay Kuramoto model to address this question.
Utilizing a recently reported analytical approach to the Kuramoto dynamics
[4], we introduce a complex-valued delay operator. This operator shapes the
dynamics of the Kuramoto system into waves traveling across the network. The
combination of this delay operator and the adjacency matrix determines these
dynamics through their effect on eigenvalues in the complex plane, thus
providing mechanistic insights into the effect of heterogeneous time delays.
The approach introduced here offers a mathematical description for the
dynamics of time-delayed networks, an important open problem in physics [5]
with many applications in neuroscience [6], engineering [7], and technology
[8]. In general, approaches to systems with heterogeneous time delays center
on numerical simulations, and no coherent analytical approach currently exists
[9, 10]. Importantly, while this question first arose from observations of
neural dynamics in the human cortex during sleep, the delay operator we
introduce here is general to studying the effect of distributed time delays in
networks at finite scales, potentially allowing insight into these dynamics in
a broad range of systems [11, 12, 13].

##  II Delay operator

We start with the standard Kuramoto model (KM) [14, 15, 16] and then consider
the model with distance-dependent time delays [17, 18, 19]. The original KM on
a general network of $N$ nodes is defined by:

$\dot{\theta}_{i}(t)=\omega_{i}+\epsilon\sum_{j=1}^{N}A_{ij}\sin(\theta_{j}(t)-\theta_{i}(t))\,,$
(1)

where $\theta_{i}\in[-\pi,\pi)$ represents the state variable (phase) of
oscillator $i$ at time $t$, $\omega_{i}$ is the intrinsic angular frequency,
$\epsilon$ scales the coupling strength, and $A_{ij}\in\\{0,1\\}$ represents
the elements of the adjacency matrix. The coupling of two connected
oscillators $i$ and $j$ causes their phases to attract [15, 16, 20, 21].

Time delays have been observed to be an important mechanism underlying the
generation of traveling waves in the brain [22, 23, 13, 24]. With this in
mind, we consider a time-delay Kuramoto model (dKM) with delays $\tau_{ij}$
that depend on the distance between two oscillators $i$ and $j$:

$\dot{\theta}_{i}(t)=\omega_{i}+\epsilon\sum_{j=1}^{N}A_{ij}\sin\Big{(}\theta_{j}(t-\tau_{ij})-\theta_{i}(t)\Big{)}\,.$
(2)

The delay operator approach we introduce here generalizes to arbitrary
adjacency matrices. In order to demonstrate this approach, we start by
considering an undirected ring graph $\mathfrak{G}_{RG}$, where $N=100$ nodes
are arranged on a one-dimensional ring with periodic boundary conditions. Each
node in $\mathfrak{G}_{RG}$ is connected to the $k=25$ nearest neighbors in
each direction, and $A_{ij}\in\\{0,1\\}$ is $1$ if oscillators $i$ and $j$ are
connected, and $0$ otherwise. The time delay $\tau_{ij}=d_{ij}/\nu$ between
two nodes $i$ and $j$ grows linearly with distance ($d_{ij}$) with respect to
the periodic boundary conditions on the ring ($d_{ij}=\min(|i-j|,N-|i-j|)$).
For the parameters chosen in this work, the time delays range from
approximately $2$ to $62$ ms, a timescale relevant to neural dynamics [25, 26,
10]. We consider the case where all oscillators have the same frequency of
$10$ Hz ($\omega=20\pi$); however, our approach can be applied to the case of
non-identical natural frequencies [27].

The time-delay term $\theta_{j}(t-\tau_{ij})$ can be approximated by
$\theta_{j}(t)-\omega\tau_{ij}$ [18, 17, 24]. Using this approximation, in
combination with the algebraic approach to the Kuramoto dynamics [4, 28], we
introduce a delay operator, which provides analytical insight into how
heterogeneous time delays can create specific, sophisticated spatiotemporal
structures in the resulting nonlinear dynamics. Applying this approximation to
Eq. (2), we arrive at an equation that captures the time-delay dynamics in the
dKM in heterogeneous phase lags [18, 17, 24]. We can then use our algebraic
approach to the Kuramoto dynamics and arrive at (see Appendix – Sec. V.1 – for
details)

$\bm{x}(t)=e^{\mathrm{i}\omega t}e^{t\bm{W}}\bm{x}(0)\,,$ (3)

where $\bm{x}\in\mathbb{C}^{N}$, and the matrix $\bm{W}$ is given by

$\bm{W}=\epsilon e^{-\mathrm{i}\bm{\eta}}\circ\bm{A}\,,$ (4)

where $\circ$ represents the Hadamard (elementwise) product. This matrix has
information about the coupling strength $\epsilon$, the time delays
$\bm{\eta}=\omega\bm{\tau}$ present in the original dKM, and the connection
scheme of the system on $\bm{A}$. In previous work, we have shown this
complex-valued equation, when evaluated through the procedure described below,
precisely captures the trajectories of the original, nonlinear Kuramoto model
[28]. We now show that this approach generalizes to the case of heterogeneous
time delays.

With this approach, we have two dynamical systems: the original, nonlinear KM
and a complex-valued system with the explicit solution in Eq. (3) (details on
the derivation can be found in [28] and in the Appendix - Sec. V.1). In the
complex-valued system, $\bm{x}\in\mathbb{C}^{N}$ has elements
$x_{i}(t)\in\mathbb{C}$ whose argument we compare with the numerical solution
of the original Kuramoto model with heterogeneous time delays (dKM)
$\theta_{i}(t)\in\mathbb{R}$ (obtained by Euler integration of Eq. (2) with
high temporal precision). That is, $\mathrm{Arg}[x_{i}(t)]$ is compared with
$\theta_{i}(t)$. When initialized with unit-modulus initial conditions
$|x_{i}(0)|=1\leavevmode\nobreak\ \text{for all}\leavevmode\nobreak\ i$, with
arguments $\mathrm{Arg}[x_{i}(0)]$ that match the initial phases
$\theta_{i}(0)$ in the original dKM, the trajectories in the original and
complex-valued KM correspond for a non-trivial window of time [4]. As
mentioned above, in [28] we found that iterating the explicit expression Eq.
(3) in a specific manner produces trajectories in the complex-valued system
that precisely match those in the original, nonlinear Kuramoto model.
Specifically, we can evaluate:

$\bm{x}(t+\varsigma)=\Lambda\big{[}e^{\mathrm{i}\omega\varsigma}e^{\varsigma\bm{W}}\bm{x}(t)\big{]}\,,$
(5)

where $\varsigma$ is small but finite,
$t\in[0,\varsigma,2\varsigma,\cdots,n\varsigma]$, and $\Lambda$ represents an
elementwise operator mapping the modulus of each state vector element
$x_{i}(t)$ to unity. This approach represents an iterative analytical
procedure, defined by the application of the linear matrix exponential and
$\Lambda$. Note that Eq. (5) propagates the solution at discrete time
intervals defined by $\varsigma$, Eq. (3) can be applied within intervals
defined by $\varsigma$, and $\varsigma>dt$. Critically, while this iterative
procedure does not represent a closed-form, all-time solution for the dynamics
of the original nonlinear Kuramoto system, all evolution of the arguments
$\mathrm{Arg}[x_{i}]$ (which, again, correspond with
$\theta_{i}(t)\leavevmode\nobreak\ \forall\leavevmode\nobreak\ i$ in the
original KM) is governed under the linear matrix exponential operator, and it
is clear that the elementwise $\Lambda$ operator only changes the moduli. In
this work, we show that this approach applies also in the case of
heterogeneous time delays and provides analytical insight into how distance-
dependent time delays create specific spatiotemporal patterns.

##  III Results

We first study phase synchronization in networks with (dKM) and without (KM)
time-delays on $\mathfrak{G}_{RG}$, as a function of the coupling strength
$\epsilon$ (Fig. 1). We use the Kuramoto order parameter:

$R(t)=\frac{1}{N}\left|\sum_{j=1}^{N}e^{\mathrm{i}\theta_{j}(t)}\right|,$ (6)

and its time average $\langle R\rangle$ for 10-second simulations to measure
the level of phase synchronization. As the coupling strength $\epsilon$
increases in the non-delayed case (original KM and complex-valued KM),
$\langle R\rangle$ begins at a low value and increases until approaching unity
(representing phase synchronization).

Figure 1: Synchronization level for non-delayed and delayed networks. The
time-average Kuramoto order parameter $\langle R\rangle$ is plotted as a
function of the coupling strength $\epsilon$ for the non-delayed case (blue
dots: original KM, orange diamonds: complex-valued KM) and for the delayed
case (red squares: original dKM, green triangles: complex-valued dKM). Each
dot represents one 10-second simulation with random initial conditions
($\mathcal{U}(-\pi,\pi)$), which are the same for the complex-valued case and
for the numerical simulation at each point. Figure 2: Analytical and geometric
view on the effect of time delays. The spatiotemporal dynamics of the system
is represented in color-code, where the phase of each oscillator is plotted as
a function of time for (a) the original KM, (b) the original dKM, (c) and the
complex-valued dKM. Dark colors represent phases close to $-\pi$ and light
colors phases close to $\pi$. (a) Without delay, the network transitions to
phase synchronization, which is represented by the horizontal lines. The
effect of the delay, however, induces wave patterns in the system, whose
dynamics are represented (b) in the original dKM and also captured (c) by the
complex-valued model. (d) These dynamical characteristics are corroborated by
the Kuramoto order parameter $R(t)$. (e) The eigenmodes offer a geometric
perspective to such dynamics, where the waves are represented by a single
eigenmode contribution ($3^{\mathrm{rd}}$ mode in this case). (f) The
eigenvalues of $\bm{W}$ (delayed) and $\epsilon\bm{A}$ (non-delayed) provide
further analytical insights into the effect of the delay in the system: it
rotates the eigenvalues in the complex plane, which allows the system to
access different modes. In the non-delayed case, the leading eigenvalue (in
real part) is associated with an eigenvector with a zero phase difference
configuration ($1^{\mathrm{st}}$ mode). In the delayed case, otherwise, there
are two leading eigenvalues that are associated to the eigenvectors
$\bm{v}_{3}$ and $\bm{v}_{99}$, which have phase configurations representing
traveling waves.

In the case with heterogeneous time delay (original dKM and complex-valued
dKM), the order parameter remains low (red squares and green triangles),
reflecting the fact that time delays induce a range of spatiotemporal
patterns, as observed previously [29, 17, 18, 30, 31, 32, 33]. Here, we
observe that the complex-valued model is able to capture the average dynamics
that the original Kuramoto model depicts, for both the non-delayed and delayed
cases, for different coupling strength across different initial conditions
(Fig. 1).

We next study dynamics in the KM and dKM considering an individual
realization, for a fixed coupling strength ($\epsilon=0.5$), and compare the
dynamics of the original dKM to the evaluation of the complex-valued approach.
Without time delays, the original KM exhibits a quick transition from random
initial conditions to a phase synchronized state (horizontal lines, Fig. 2a).
With time delays, however, phase synchronization is not reached, and the
original dKM exhibits a transition from random initial conditions to a
traveling wave state (diagonal structures, Fig. 2b). The evaluation of the
complex-valued dKM captures both the transient dynamics and traveling wave
state exhibited in the original dKM (Fig. 2c), as well as the dynamics of the
Kuramoto order parameter $R(t)$ (Fig. 2d).

Our approach to systems with heterogeneous time delay provides insight into
the mechanism for these dynamics in terms of the spectrum of $\bm{W}$ – Eq.
(4). If $\bm{A}$ and $\bm{\tau}$ are circulant, $\bm{W}$ is also circulant
(see Appendix – Sec. V.2), hence $\bm{W}$ and $\bm{A}$ share the same
eigenvectors (which form an orthonormal basis). We can then write Eq. (3)
using the eigenspectrum of $\bm{W}$, which results in
$\bm{x}(t)=e^{\mathrm{i}\omega
t}(\alpha_{1}e^{\lambda_{1}t}\bm{v}_{1}+\cdots+\alpha_{N}e^{\lambda_{N}t}\bm{v}_{N})$,
where $\alpha_{i}$ can be written in terms of initial conditions. Importantly,
we can also write Eq. (5) in a similar fashion, which results in
$\bm{x}(t+\varsigma)=\Lambda\big{[}e^{\mathrm{i}\omega\varsigma}(\alpha_{1}e^{\lambda_{1}\varsigma}\bm{v}_{1}+\cdots+\alpha_{N}e^{\lambda_{N}\varsigma}\bm{v}_{N})\big{]}$,
where $\alpha_{i}$ can again be written in terms of the state of the system at
time $t\in[0,\varsigma,2\varsigma,\cdots,n\varsigma]$. Thus, while it is in
general a very difficult problem to understand the dynamics of nonlinear
networks in terms of eigenspectra, this approach provides a unique insight
into the connection between the spectrum of $\bm{W}$ – Eq. (4) – and the
spatiotemporal dynamics of the nonlinear oscillator network – Eq. (2).
Critically, our approach uses familiar mathematical techniques from linear
algebra matrix theory in a distinct way: while previous approaches in
nonlinear dynamics have sought to describe the dynamics using the spectrum of
the Laplacian matrix [34, 35, 36], the focus on the complex-valued system in
our approach enables the insight that the argument of the eigenvectors of the
matrix $\bm{W}$ provides analytical predictions about the resulting nonlinear
dynamics.

Following this idea, Fig. 2e shows the eigenmode contributions, here
represented by $\log{|\mu_{i}|}$, as a function of time, for the dynamics in
Fig. 2c. Here, the eigenmode contributions are given by the projection of the
complex-valued approach solution $\bm{x}(t)$ onto the eigenvectors of
$\bm{W}$. The eigenmode contributions are obtained as
$\mu_{k}(t)=\langle\bm{x}(t),\bm{v}_{k}\rangle$, where $\langle.\rangle$
denotes the standard complex inner product. Figure 2e shows that, when the
network exhibits incoherent dynamics, the eigenmode contributions remain
uniform across $\mu_{i}$. When the traveling wave pattern is reached, on the
other hand, the $3^{\mathrm{rd}}$ eigenmode becomes dominant (note the log
scale). These results demonstrate that the change from incoherent dynamics to
a traveling wave can be understood quite directly through the geometry of the
eigenmodes. Further, in the case of circulant networks, we can evaluate
eigenvalues and eigenvectors analytically using the circulant diagonalization
theorem (CDT) [37], and in this case, the $1^{\mathrm{st}}$ eigenvector
represents the solution where all oscillators have the same phase (phase
synchronization), and higher modes represent wave patterns, given by Fourier
modes (see Appendix - Sec. V.3).

The effect of heterogeneous time delays on the dynamics of the dKM can be
understood through the geometry of eigenvalues in the complex plane. Figure 2f
illustrates the eigenvalues of $\epsilon\bm{A}$ (non-delayed) and $\bm{W}$
(delayed). While the non-delayed case (blue line and dots) has purely real
eigenvalues, the effect of the heterogeneous time delays (red line and
squares) can be understood in our framework in terms of the Hadamard
(elementwise) product of the delay operator $\bm{\tau}$ and $\bm{A}$ (see Eq.
(4), and Appendix - Sec. V.2). The effect of this operation is to provide a
specific rotation of the eigenvalues in the complex plane. This rotation
allows the system to access higher modes and, therefore, to exhibit different
traveling wave patterns. Further, the rotation is not the same for all
eigenvalues because the delays are heterogeneous. In this particular case, the
rotation leads to eigenvalues associated to the $3^{\mathrm{rd}}$ and
$99^{\mathrm{th}}$ modes to have the largest real part, allowing the system to
reach traveling wave states associated with the $3^{\mathrm{rd}}$ and
$99^{\mathrm{th}}$ modes. In the particular example of Fig. 2, the network
evolves to a wave given by the $3^{\mathrm{rd}}$ mode, but different (random)
initial conditions can either evolve to the dynamics described by the
$3^{\mathrm{rd}}$ or $99^{\mathrm{th}}$ mode [38]. Moreover, when different
time delays are considered, different modes can be dominant, and therefore the
system evolves to a different wave pattern [39].

Figure 3: Analytical predictions of specific wave patterns. (a) The phase
configuration for the original dKM (blue line) matches the argument
(elementwise) of the $3^{\mathrm{rd}}$ eigenvector $\mathrm{Arg}[\bm{v}_{3}]$
\- analytical prediction - (red dotted line). A representation on the circle
using a color-code reveals the wave pattern (right). (b) Different initial
conditions lead to the wave pattern that matches the argument of the
$99^{\mathrm{th}}$ eigenvector $\mathrm{Arg}[\bm{v}_{99}]$. These waves can
propagate either counterclockwise (negative) or clockwise (positive). (c) With
random initial conditions, due to the dominance of two eigenvalues
($3^{\mathrm{rd}}$ and $99^{\mathrm{th}}$), the system exhibits waves
propagating in both directions - with approximately half of the initial
conditions evolving to each direction (top right). (d) With biased initial
conditions, starting from ${\mathrm{Arg}[\bm{v}_{3}]}$ (red line, bottom
right) and adding uniform random phases $0.8(\mathcal{U}(-\pi,\pi))$, we
obtain a preferred direction of propagation.

We can now uncover how the combination of network structure, time delays, and
node state can create specific spatiotemporal patterns. By using our delay
operator approach, we can analytically predict the specific pattern to which
the original dKM evolves. Figure 3a shows the wave pattern given by $\theta$
obtained from the original dKM (blue line) and the argument (elementwise) of
the $3^{\mathrm{rd}}$ eigenvector (red line), which predicts the observed
dynamics [40]. In this case, phases increase in the clockwise direction around
the ring, which we define to be the positive direction ($+1$). It is important
to note that, in our approach, the argument of each eigenvector element
$(\mathrm{Arg}[(\bm{v}_{k})_{i}]\,\forall i\in[1,N])$ directly relates with
the phase offset in the resulting network dynamics. Because of the
correspondence between trajectories in the complex-valued model and the
original dKM, this approach creates a direct link between eigenvectors of the
adjacency matrix and the specific spatiotemporal dynamics that result. For the
dynamics in Fig. 3a, the eigenmode contribution is given by $\mu_{3}$ (see
Fig. 2e), and the phase configuration matches the argument of $\bm{v}_{3}$. In
the example considered here, two eigenvalues are dominant (i.e. having the
largest real part) - $\lambda_{3}$ and $\lambda_{99}$ (Fig. 2f). Different
initial conditions can thus evolve to the phase pattern given by the
$99^{\mathrm{th}}$ mode, which is predicted by $\bm{v}_{99}$ (Fig. 3b). In
this case, the spatial frequency is the same as observed in the previous case,
but the direction of the wave pattern is the opposite [41]. These results show
a clear connection between the spectrum of the network (described by $\bm{W}$)
and the dynamics on the original dKM, where the wave pattern (solution) can be
described by the phase configuration of the eigenvector associated to the
dominant mode.

We take counterclockwise increases in phase to be in the negative direction,
and clockwise increases to be positive. Because the network considered here
has two dominant eigenvalues equal in their real parts, random initial
conditions evolve equally either to the phase pattern of $\bm{v}_{3}$ or
$\bm{v}_{99}$ in individual simulations (Fig. 3c). To quantify the
spatiotemporal dynamics, the spatial frequency, and the direction of
propagation, we compare the phases obtained from the original dKM and the
argument of the eigenvectors of $\bm{W}$. Specifically, we evaluate:

$\rho^{(k)}(t)=\left|\frac{1}{N}\sum_{i}^{N}e^{\mathrm{i}\theta_{i}(t)}e^{-\mathrm{i}\mathrm{Arg}[(\bm{v}_{k})_{i}]}\right|,$
(7)

where $\theta_{i}(t)$ is the phase of the oscillators $i$ at time $t$ obtained
from the original dKM, $N$ is the number of oscillators in the network,
$\mathrm{i}$ is the imaginary unit, and $\bm{v}_{k}$ is the $k^{\mathrm{th}}$
eigenvector of $\bm{W}$. Here, we use $\bm{v}_{3}$, and $\rho^{(k)}=1$ means
the phase configuration of the network given by the $\bm{\theta}(t)$ is the
same as the one given by the argument of the eigenvector $\bm{v}_{k}$. In the
case shown in Fig. 3c, approximately half of the simulations evolve to the
positive direction, indicating the dynamics matches the argument of
$\bm{v}_{3}$, and approximately half evolve to the negative, indicating the
dynamics is given by the argument of $\bm{v}_{99}$. A small fraction of
initial conditions exhibit inner products approximately $\pm 0.5$,
corresponding to a wave with a different spatial frequency.

Using the insights from this approach, we can now design initial conditions
that generate waves in a preferred direction. To do this, we started from the
phase pattern specified by $\bm{v}_{3}$ and randomized the phases by nearly a
full cycle ($0.8\leavevmode\nobreak\ \mathcal{U}[-\pi,\pi]$, then wrapped in
$[-\pi,\,\pi]$). While this initial condition is nearly random (Fig. 3d,
bottom right, where the red line represents ${\mathrm{Arg}[\bm{v}_{3}]}$;
compare with Fig. 3c, bottom right), nearly all simulations evolve to the
positive direction. These results demonstrate that the combination of
connectivity, time delays, and network state can generate specific
spatiotemporal patterns in oscillator networks – here, traveling waves with a
chirality in a preferred direction.

The framework for systems with heterogeneous time delays introduced in the
work generalizes to many types of networks. This approach can be applied to
very sophisticated networks obtained from experimental data. In particular,
this approach can successfully predict traveling wave patterns arising in an
oscillator network based on connectivity in the human brain. Figure 4
illustrates simulations and the analytical prediction resulting from our
approach for networks where the connectivity data is based on the Human
Connectome Project (HCP) [42]. In this case, $N=998$ cortical regions are
given at a point in 3-space, with connections between areas derived from
neuroimaging data. Connection weights between regions are determined by the
number of fibers [42, 43], which we use to build the adjacency matrix
$\bm{A}$. Here, the coupling strength is scaled with $\epsilon=200$, and the
initial conditions for each analysis are given by random phases $[-\pi,\pi]$.
Further, time delays are obtained by $\tau_{ij}=d_{ij}/\nu$, where the
distances $d_{ij}$ are determined by the average length of these fibers, and
the known axonal conduction speed is given by $\nu=5$ m/s [44]. The dynamics
of each node is represented by the Kuramoto model, given either by Eq. (1) in
the nondelayed case and by Eq. (2) in the delayed case. The natural frequency
of each oscillator is given by $10$ Hz (simulating, for example, a specific
drive from the thalamus). Using the delay operator, we construct the matrix
$\bm{W}$ for these systems – Eq. (4) – which allows us to obtain analytical
predictions of the spatiotemporal patterns that emerge.

Figure 4: Analytical predictions of spatiotemporal patterns in brain networks.
We use our approach to investigate networks based on the Human Connectome
Project (HCP) [42]. We consider the case without delay in the coupling between
nodes and also the case with distance dependent delays (heterogeneous delay).
We use our delay operator to create the matrix $\bm{W}$, which allows
analytical predictions of the dynamics. We show the phase of each node, given
by the Kuramoto model, using a color-code (dark colors are values close to
$-\pi$, and light colors are values close to $\pi$). (a) In the case without
delay, the argument of the leading eigenvector depicts zero phase difference,
which predicts phase synchronization. (b) We then study the numerical
simulation for the network without delay, given by Eq. (1), which shows a
phase synchronized state. (c) In the case with heterogeneous time delays, the
argument of the leading eigenvector shows a phase offset from bottom left to
the top right (in the projection), which predicts a wave pattern. (d) We
perform the numerical simulation with the delayed Kuramoto model, given by Eq.
(2), and the network depicts the wave pattern that is predicted by our
approach. This shows that we are able to predict the dynamics observed in the
simulations using our delay operator.

First, we consider the case without time delays, where $\tau_{ij}=0$. We then
obtain the eigenspectrum of the matrix $\bm{W}$ and plot the argument
(elementwise) of the eigenvector associated with the leading eigenvalue (Fig.
4a). In this case, this eigenvector shows a zero phase difference across
nodes, predicting phase synchronization. We then perform the numerical
simulation of the Kuramoto model (without delay), given by Eq. (1) and plot
the phase of each node in color-code (Fig. 4b), where we observe a phase
synchronized behavior [45].

On the other hand, when we consider time delays in the interaction between
cortical areas, the scenario is different. In this case, the argument
(elementwise) of the eigenvector associated with the leading eigenvalue
depicts a phase offset increasing from the bottom left to top right (in this
projection), predicting a wave propagating along that direction (Fig. 4c). We
then perform numerical simulations of the Kuramoto model with heterogeneous
time delays – Eq. (2) and we observe the wave pattern that is predicted by our
approach, as shown in Fig. 4d [46].

This example now clearly demonstrates the advantage of this analytical
approach: when we numerically evaluate the eigenspectrum of $\bm{W}$ in this
case, the leading eigenvector for the case without delays predicts phase
synchrony, while the leading eigenvector for the case with delays predicts the
precise wave pattern observed in the simulation. This result shows that our
approach is able to predict the spatiotemporal pattern that results from
connectivity and time delays in a highly relevant, real-world case.

##  IV Conclusion

In this work, we have introduced an analytical approach to the dynamics of
nonlinear oscillator networks with heterogeneous time delays, an important
open problem in physics with many potential applications. The advance in this
work is based on an algebraic approach to the Kuramoto model introduced in
[28]. Importantly, the flexibility of this framework allowed us to introduce a
delay operator, which provides rigorous analytical predictions for the
specific traveling wave patterns induced by distance-dependent time delays.
Using this approach, we can explain the effect of time delays in terms of a
rotation of the eigenvalues of the matrix describing the system, which
provides a clear and precise way to understand heterogeneous time delays in
terms of the geometry of eigenmodes. Our approach therefore allows analytical
predictions for the specific spatiotemporal patterns exhibited by the original
dKM.

This framework allows us to understand how the combination of isotropic
connectivity and time delays can produce traveling waves propagating in a
preferred direction, as observed in experimental data [1]. Importantly, while
this question first arose in our study of neural dynamics in human cortex
during sleep, the approach we have introduced here is general to networks of
oscillators at finite scales. The results shown in this work, together with
the results in [28, 4], represent a coherent and general framework for
nonlinear oscillator networks.

The central advance of this framework is to consider the dynamics in an
individual simulation, taking into account both the initial conditions and the
specific connectivity pattern in the network. This framework thus provides an
opportunity to connect an individual adjacency matrix, for example a single
network taken from experimental data or a single realization of a random graph
model, to the specific spatiotemporal pattern that results in a simulation.
This approach has important potential applications, for example in linking an
experimentally reconstructed brain network to dynamics and computation in a
neural system [47] or in linking the connections in a large-scale power grid
to potential large and transient disruptions [48, 49]. In this work, we have
generalized this framework to systems with heterogeneous time delays, which
demonstrates the utility of this algebraic, operator-based approach to
nonlinear dynamical systems at finite scales.

###### Acknowledgements.

This work was supported by BrainsCAN at Western University through the Canada
First Research Excellence Fund (CFREF), the NSF through a NeuroNex award
(#2015276), the Natural Sciences and Engineering Research Council of Canada
(NSERC) grant R0370A01, ONR N00014-16-1-2829, NIH EB009282, NIH EB026899, the
Swartz Foundation, SPIRITS 2020 of Kyoto University, Compute Ontario
(computeontario.ca), Compute Canada (computecanada.ca), and the Western
Academy for Advanced Research. J.M. gratefully acknowledges the Western
University Faculty of Science Distinguished Professorship in 2020-2021. R.C.B
gratefully acknowledges the Western Institute for Neuroscience Clinical
Research Postdoctoral Fellowship. G.B.B gratefully acknowledges the Canadian
Open Neuroscience Platform (Graduate Scholarship), the Vector Institute
(Postgraduate Affiliate), and NSERC (Canada Graduate Scholarship - Doctoral).

## Code availability

An open-source code repository for this work is available on GitHub:
http://mullerlab.github.io.

##  V Appendix

###  V.1 The complex-valued approach

We consider the Kuramoto model with heterogeneous time delays described by Eq.
(2) and then use the approximation given by
$\theta_{j}(t-\tau_{ij})\approx\theta_{j}(t)-\omega\tau_{ij}$ [18, 17, 24],
which leads to

$\dot{\theta}_{i}(t)=\omega+\epsilon\sum_{j=1}^{N}A_{ij}\sin\Big{(}\theta_{j}(t)-\theta_{i}(t)-\eta_{ij}\Big{)},$
(8)

where $\eta_{ij}=\omega\tau_{ij}$.

Based on [4, 28], we introduce the complex-valued approach to the Kuramoto
model described by Eq. (8). To do that, we introduce a new dynamical system,
described by the variable $\psi\in\mathbb{C}$:

$\begin{split}\dot{\psi}_{i}(t)=\omega+\epsilon\sum_{j=1}^{N}A_{ij}\Big{[}\sin\Big{(}\psi_{j}(t)-\psi_{i}(t)-\eta_{ij}\Big{)}\\\
-\mathrm{i}\cos\Big{(}\psi_{j}(t)-\psi_{i}(t)-\eta_{ij}\Big{)}\Big{]}\,.\end{split}$
(9)

Next, multiplying both sides by $\mathrm{i}$ and applying Euler’s formula
yields

$\mathrm{i}\dot{\psi}_{i}(t)=\mathrm{i}\omega+\epsilon
e^{-\mathrm{i}\psi_{i}(t)}\sum_{j=1}^{N}A_{ij}e^{\mathrm{i}\psi_{j}(t)}e^{-\mathrm{i}\eta_{ij}}.$
(10)

We define $\bm{W}$ as:

$\bm{W}=\epsilon e^{-\mathrm{i}\bm{\eta}}\circ\bm{A}\,,$ (11)

where $\circ$ represents the Hadamard product (or element-wise product), and
$\eta_{ij}=\omega\tau_{ij}$. This results in the following matrix form of Eq.
(10):

$\bm{\dot{\psi}}(t)=\bm{\omega}+\frac{1}{\mathrm{i}}\,\text{diag}[e^{-\mathrm{i}\bm{\psi}(t)}]\,\bm{W}e^{\mathrm{i}\bm{\psi}(t)}\,,$
(12)

where we note explicitly that $\bm{\psi}=[\psi_{1},\cdots,\psi_{N}]^{T}$,
$\bm{\dot{\psi}}=[\dot{\psi}_{1},\cdots,\dot{\psi}_{N}]^{T}$, and
$\bm{\omega}=[\omega,\cdots,\omega]^{T}$. Furthermore, we can write the
previous equation as:

$\frac{d}{dt}e^{\mathrm{i}\bm{\psi}(t)}=\bigg{(}\text{diag}[\mathrm{i}\bm{\omega}]+\bm{W}\bigg{)}e^{\mathrm{i}\bm{\psi}(t)}\,.$
(13)

Lastly, letting $\bm{x}(t)=e^{\mathrm{i}\bm{\psi}(t)}$, we have

$\bm{\dot{x}}(t)=\bigg{(}\text{diag}[\mathrm{i}\bm{\omega}]+\bm{W}\bigg{)}\bm{x}(t)\,,$
(14)

whose general solution is

$\bm{x}(t)=e^{\mathrm{i}\omega t}e^{t\bm{W}}\bm{x}(0)\,.$ (15)

In this work, the dynamics of the complex-valued approach is studied by
considering the elementwise argument of $\bm{x}(t)$, i.e.
$\mathrm{Arg}[x_{i}(t)]\leavevmode\nobreak\ \forall\leavevmode\nobreak\
i\in[1,N]$. As shown in [28], when $\frac{|x_{j}|}{|x_{i}|}\approx 1$, the
dynamics of $\mathrm{Arg}[\bm{x}(t)]$ precisely matches the trajectories of
the Kuramoto model given by Eq. (8). This allows us to use the eigenspectrum
of $\bm{W}$ to understand and predict the dynamics of the Kuramoto model with
heterogeneous time delays.

###  V.2 Circulant networks and Hadamard product

The definition of the Hadamard product can be described as follows:

######  Definition 1.

Let $\bm{A},\bm{B}$ be two $n\times n$ matrices. The Hadamard product
$\bm{A}\circ\bm{B}$ is a matrix of dimension $n\times n$ with elements given
by

$(\bm{A}\circ\bm{B})_{ij}=(\bm{A})_{ij}(\bm{B})_{ij}.$

For a complex number $\lambda$, we also define $e^{\circ(\lambda\bm{A})}$ to
be the matrix of dimension $n\times n$ with elements given by

$(e^{\circ\lambda\bm{A}})_{ij}=e^{\lambda\bm{A}_{ij}}.$

We have the following observation.

######  Proposition 1.

Let $\bm{A},\bm{B}$ be two circulant matrices. Then

  1. 1.

$\bm{A}\circ\bm{B}$ is a circulant matrix.

  2. 2.

$e^{\lambda\circ\bm{A}}$ is a circulant matrix.

###### Proof.

Assume that $\bm{A}=\text{circ}(\bm{a}),\bm{B}=\text{circ}(\bm{b})$ with
$\bm{a}=(a_{1},a_{2},\ldots,a_{n})$ and $\bm{b}=(b_{1},b_{2},\ldots,b_{n}).$
Then we can see that

$\bm{A}\circ\bm{B}=\text{circ}((a_{1}b_{1},a_{2}b_{2},\ldots,a_{n}b_{n})),$

and

$e^{\circ\lambda\bm{A}}=\text{circ}((e^{\lambda a_{1}},e^{\lambda
a_{2}},\ldots,e^{\lambda a_{n}})).$

Therefore, we conclude that both $\bm{A}\circ\bm{B}$ and
$e^{\lambda\circ\bm{A}}$ are circulant. ∎

###  V.3 The Circulant Diagonalization Theorem

In the case of circulant networks, we can use the circulant diagonalization
theorem (CDT) to obtain the eigenspectrum of the adjacency matrix analytically
[37]. In our work, both the non-delayed network $\epsilon\bm{A}$ and the
delayed one $\bm{W}$ are circulant (see Prop. 1).

Figure 5: A graphical representation of the matrix with the phase
configuration of the eigenvectors of $\bm{W}$. The $k^{\mathrm{th}}$ column is
the color-coded argument (elementwise) of the $k^{\mathrm{th}}$ eigenvector.

The CDT states that all circulant matrices, say $\bm{H}=\text{circ}(\bm{h})$,
where $\text{circ}(\bm{h})$ is the circulant matrix constructed from the
generating vector $\bm{h}=(h_{1},\cdots,h_{N})$, are diagonalized by the same
unitary matrix $\bm{U}$ with components

$U_{ks}=\frac{1}{\sqrt{N}}\exp\left[-\frac{2\pi\mathrm{i}}{N}(k-1)(s-1)\right],$
(16)

where $k,s\in[1,N]$ and that the $N$ eigenvalues are given by

$E_{k}(\bm{H})=\sum\limits_{j=1}^{N}h_{j}\exp\left[-\frac{2\pi\mathrm{i}}{N}(k-1)(j-1)\right]\,.$
(17)

We let Eq. (17) determine the ordering of the eigenvalues throughout this
work. The argument of the eigenvectors associated with these eigenvalues
correspond to the columns of the discrete Fourier transform (DFT) matrix,
which range from low to high spatial frequencies.

Figure 5 shows the argument of the eigenvectors in color-code. Here,
$\mathrm{Arg}[(\bm{v}_{1})_{i}]=0\leavevmode\nobreak\
\forall\leavevmode\nobreak\ i\in[1,N]$ (as shown in Fig. 2), which represents
zero phase difference across oscillators, or phase synchronization. The other
eigenvectors represent Fourier modes (waves) with different spatial
frequencies. Figure 2 shows the cases of the eigenvectors $\bm{v}_{3}$ and
$\bm{v}_{99}$.

## References

  * Muller  _et al._ [2016] L. Muller, G. Piantoni, D. Koller, S. S. Cash, E. Halgren, and T. J. Sejnowski, Rotating waves during human sleep spindles organize global patterns of activity that repeat precisely through the night, Elife 5, e17267 (2016).
  * Sejnowski and Destexhe [2000] T. J. Sejnowski and A. Destexhe, Why do we sleep?, Brain research 886, 208 (2000).
  * Girard  _et al._ [2001] P. Girard, J. M. Hupé, and J. Bullier, Feedforward and feedback connections between areas v1 and v2 of the monkey have similar rapid conduction velocities, Journal of Neurophysiology 85, 1328 (2001).
  * Muller  _et al._ [2021] L. Muller, J. Mináč, and T. T. Nguyen, Algebraic approach to the kuramoto model, Physical Review E 104, L022201 (2021).
  * Atay [2010] F. M. Atay,  _Complex time-delay systems: theory and applications_ (Springer, 2010).
  * Schröter  _et al._ [2017] M. Schröter, O. Paulsen, and E. T. Bullmore, Micro-connectomics: probing the organization of neuronal networks at the cellular scale, Nature Reviews Neuroscience 18, 131 (2017).
  * Papachristodoulou  _et al._ [2010] A. Papachristodoulou, A. Jadbabaie, and U. Münz, Effects of delay in multi-agent consensus and oscillator synchronization, IEEE transactions on automatic control 55, 1471 (2010).
  * Liao and Wang [2000] T. Liao and F. Wang, Global stability for cellular neural networks with time delay, IEEE Transactions on neural networks 11, 1481 (2000).
  * Tewarie  _et al._ [2019] P. Tewarie, R. Abeysuriya, A. Byrne, G. C. O’Neill, S. N. Sotiropoulos, M. J. Brookes, and S. Coombes, How do spatially distinct frequency specific meg networks emerge from one underlying structural connectome? the role of the structural eigenmodes, NeuroImage 186, 211 (2019).
  * Petkoski and Jirsa [2022] S. Petkoski and V. K. Jirsa, Normalizing the brain connectome for communication through synchronization, Network Neuroscience , 1 (2022).
  * Lee  _et al._ [2009] W. S. Lee, E. Ott, and T. M. Antonsen, Large coupled oscillator systems with heterogeneous interaction delays, Physical Review Letters 103, 044101 (2009).
  * Papachristodoulou and Jadbabaie [2006] A. Papachristodoulou and A. Jadbabaie, Synchonization in oscillator networks with heterogeneous delays, switching topologies and nonlinear dynamics, in  _Proceedings of the 45th IEEE Conference on Decision and Control_ (IEEE, 2006) pp. 4307–4312.
  * Roberts  _et al._ [2019] J. A. Roberts, L. L. Gollo, R. G. Abeysuriya, G. Roberts, P. B. Mitchell, M. W. Woolrich, and M. Breakspear, Metastable brain waves, Nature Communications 10, 1 (2019).
  * Kuramoto [1975] Y. Kuramoto, Self-entrainment of a population of coupled non-linear oscillators, in  _International symposium on mathematical problems in theoretical physics_ (Springer, 1975) pp. 420–422.
  * Rodrigues  _et al._ [2016] F. A. Rodrigues, T. K. D. M. Peron, P. Ji, and J. Kurths, The kuramoto model in complex networks, Physics Reports 610, 1 (2016).
  * Acebrón  _et al._ [2005] J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, The kuramoto model: A simple paradigm for synchronization phenomena, Reviews of Modern Physics 77, 137 (2005).
  * Jeong  _et al._ [2002] S. O. Jeong, T. W. Ko, and H. T. Moon, Time-delayed spatial patterns in a two-dimensional array of coupled oscillators, Physical Review Letters 89, 154104 (2002).
  * Ko and Ermentrout [2007] T. W. Ko and G. B. Ermentrout, Effects of axonal time delay on synchronization and wave formation in sparsely coupled neuronal oscillators, Physical Review E 76, 056206 (2007).
  * Petkoski  _et al._ [2016] S. Petkoski, A. Spiegler, T. Proix, P. Aram, J. Temprado, and V. K. Jirsa, Heterogeneity of time delays determines synchronization of coupled oscillators, Physical Review E 94, 012209 (2016).
  * Strogatz [2000] S. H. Strogatz, From kuramoto to crawford: exploring the onset of synchronization in populations of coupled oscillators, Physica D: Nonlinear Phenomena 143, 1 (2000).
  * Arenas  _et al._ [2008] A. Arenas, A. Diaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Synchronization in complex networks, Physics Reports 469, 93 (2008).
  * Muller  _et al._ [2018] L. Muller, F. Chavane, J. Reynolds, and T. J. Sejnowski, Cortical travelling waves: mechanisms and computational principles, Nature Reviews Neuroscience 19, 255 (2018).
  * Davis  _et al._ [2021] Z. W. Davis, G. B. Benigno, C. Fletterman, T. Desbordes, C. Steward, T. J. Sejnowski, J. Reynolds, and L. Muller, Spontaneous traveling waves naturally emerge from horizontal fiber time delays and travel through locally asynchronous-irregular states, Nature Communications 12, 1 (2021).
  * Breakspear  _et al._ [2010] M. Breakspear, S. Heitmann, and A. Daffertshofer, Generative models of cortical oscillations: neurobiological implications of the kuramoto model, Frontiers in Human Neuroscience 4, 190 (2010).
  * Cabral  _et al._ [2011] J. Cabral, E. Hugues, O. Sporns, and G. Deco, Role of local network oscillations in resting-state functional connectivity, Neuroimage 57, 130 (2011).
  * Choi and Mihalas [2019] H. Choi and S. Mihalas, Synchronization dependent on spatial structures of a mesoscopic whole-brain network, PLoS Computational Biology 15, e1006978 (2019).
  * fig [a] See Supplemental Material for Fig. S1 .
  * Budzinski  _et al._ [2022] R. C. Budzinski, T. T. Nguyen, J. Đoàn, J. Mináč, T. J. Sejnowski, and L. E. Muller, Geometry unites synchrony, chimeras, and waves in nonlinear oscillator networks, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 031104 (2022).
  * Yeung and Strogatz [1999] M. K. S. Yeung and S. H. Strogatz, Time delay in the kuramoto model of coupled oscillators, Physical Review Letters 82, 648 (1999).
  * Laing [2016] C. R. Laing, Travelling waves in arrays of delay-coupled phase oscillators, Chaos: An Interdisciplinary Journal of Nonlinear Science 26, 094802 (2016).
  * Timms and English [2014] L. Timms and L. Q. English, Synchronization in phase-coupled kuramoto oscillator networks with axonal delay and synaptic plasticity, Physical Review E 89, 032906 (2014).
  * Ermentrout and Ko [2009] B. Ermentrout and T. W. Ko, Delays and weakly coupled neuronal oscillators, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 367, 1097 (2009).
  * Ko  _et al._ [2004] T. W. Ko, S. O. Jeong, and H. T. Moon, Wave formation by time delays in randomly coupled oscillators, Physical Review E 69, 056106 (2004).
  * Grabow  _et al._ [2010] C. Grabow, S. M. Hill, S. Grosskinsky, and M. Timme, Do small worlds synchronize fastest?, EPL (Europhysics Letters) 90, 48002 (2010).
  * Sorrentino  _et al._ [2016] F. Sorrentino, L. M. Pecora, A. M. Hagerstrom, T. E. Murphy, and R. Roy, Complete characterization of the stability of cluster synchronization in complex dynamical networks, Science Advances 2, e1501737 (2016).
  * Sugitani  _et al._ [2021] Y. Sugitani, Y. Zhang, and A. E. Motter, Synchronizing chaos with imperfections, Physical Review Letters 126, 164101 (2021).
  * Davis [1979] P. Davis,  _Circulant matrices_ (John Wiley & Sons, 1979).
  * fig [b] See Supplemental Material for Fig. S2 .
  * fig [c] See Supplemental Material for Fig. S3 .
  * mov [a] See Supplemental Material for Movie S1 .
  * mov [b] See Supplemental Material for Movie S2 .
  * Hagmann  _et al._ [2008] P. Hagmann, L. Cammoun, X. Gigandet, R. Meuli, C. J. Honey, V. J. Wedeen, and O. Sporns, Mapping the structural core of human cerebral cortex, PLoS Biology 6, e159 (2008).
  * Muller  _et al._ [2014] L. Muller, A. Destexhe, and M. Rudolph-Lilith, Brain networks: small-worlds, after all?, New Journal of Physics 16, 105004 (2014).
  * Swadlow and Waxman [2012] H. A. Swadlow and S. G. Waxman, Axonal conduction delays, Scholarpedia 7, 1451 (2012).
  * mov [c] See Supplemental Material for Movie S3 .
  * mov [d] See Supplemental Material for Movie S4 .
  * Litvina  _et al._ [2019] E. Litvina, A. Adams, A. Barth, M. Bruchez, J. Carson, J. E. Chung, K. B. Dupre, L. M. Frank, K. M. Gates, K. M. Harris, _et al._ , Brain initiative: cutting-edge tools and resources for the community, Journal of Neuroscience 39, 8275 (2019).
  * Witthaut  _et al._ [2022] D. Witthaut, F. Hellmann, J. Kurths, S. Kettemann, H. Meyer-Ortmanns, and M. Timme, Collective nonlinear dynamics and self-organization in decentralized power grids, Reviews of Modern Physics 94, 015005 (2022).
  * Motter  _et al._ [2013] A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa, Spontaneous synchrony in power-grid networks, Nature Physics 9, 191 (2013).


# Revealing the learning process in reinforcement learning agents through
attention-oriented metrics

Charlotte Beylier
Max Planck Institute for Human Cognitive and Brain Sciences
Center for Scalable Data Analytics and Artificial Intelligence (ScaDS.AI)
Dresden/Leipzig, Germany
<EMAIL_ADDRESS>
&Simon M. Hofmann
Max Planck Institute for Human Cognitive and Brain Sciences
Leipzig, Germany
<EMAIL_ADDRESS>
&Nico Scherf
Max Planck Institute for Human Cognitive and Brain Sciences
Center for Scalable Data Analytics and Artificial Intelligence (ScaDS.AI)
Dresden/Leipzig, Germany
<EMAIL_ADDRESS>

###### Abstract

The learning process of a reinforcement learning (RL) agent remains poorly
understood beyond the mathematical formulation of its learning algorithm. To
address this gap, we introduce attention-oriented metrics (ATOMs) to
investigate the development of an RL agent’s attention during training. We
tested ATOMs on three variations of a Pong game, each designed to teach the
agent distinct behaviours, complemented by a behavioural assessment. Our
findings reveal that ATOMs successfully delineate the attention patterns of an
agent trained on each game variation, and that these differences in attention
patterns translate into differences in the agent’s behaviour. Through
continuous monitoring of ATOMs during training, we observed that the agent’s
attention developed in phases, and that these phases were consistent across
games. Finally, we noted that the agent’s attention to its paddle emerged
relatively late in the training and coincided with a marked increase in its
performance score. Overall, we believe that ATOMs could significantly enhance
our understanding of RL agents’ learning processes, which is essential for
improving their reliability and efficiency.

##  1 Introduction

Understanding the learning process of a deep Reinforcement Learning (RL) agent
is critical for improving its transparency and performance. However, the
mechanisms unfolding during an agent’s learning phase remain largely unknown
beyond the mathematical formulations of its algorithm and optimisation
process. Therefore, new metrics are needed to provide insights into what an RL
agent learns about its environment and its task during training.

Currently, the primary indicator of an agent’s learning progress is its
performance score derived from task-specific rewards. While necessary to
optimize an agent Watkins and Dayan [1992], Konda and Tsitsiklis [1999], Mnih
et al. [2013], this performance score provides a very limited view of an
agent’s development during learning and fails to explain or predict certain
behaviours in the testing phase. For instance, RL agents have been observed
exploiting faulty reward functions Clark and Amodei [2016] while achieving
high performance scores during training, or showing significant divergence
during the test phase despite similar training performances Cobbe et al.
[2019]. These examples underscore the limitations of this metric in fully
understanding an agent’s behaviours and development.

In this context, explainable AI (XAI) methods can complement the limitation of
the performance score by providing an insight into an agent understanding of
its environment and task through the explanation of its decisions Verma et al.
[2018], Shu et al. [2017]. However, the explanation format given by these
methods - which may be a decision tree Liu et al. [2019], a structural causal
model Madumal et al. [2020] or a set of linguistic rules Hein et al. [2017] \-
is often too complex to be tracked during learning. Methods focused on
saliency maps on the other hand, are easier to compute and are more versatile
but are mostly used in a qualitative way to explain the behaviour of trained
agent Zeiler and Fergus [2014], Zhou et al. [2016], Weitkamp et al. [2018],
Selvaraju et al. [2017]. As a results XAI techniques do not give information
during the learning phase of an agent.

In this work we introduce attention-oriented metrics (ATOMs) to gain insight
into an agent’s learning process through the development of its attention.
ATOMs is derived from saliency maps and combine the traceability of the
performance score metric with the explanatory power of XAI and quantify an
agent’s attention on the objects within its environment. Specifically, ATOMs
encompass two levels of description:

  * •

the ranked attention on individual objects (_hierarchical-attention_), and

  * •

the attention on combinations of such objects (_combinatorial-attention_).

To systematically evaluate ATOMs’ ability to give information about what an RL
agent has learned, we created three variation of a Pong game. Each variation
required the agent to learn distinct behaviour. In addition to these games, we
designed a behavioural experiment to test if the agent’s attention described
by ATOMs translated into their observed behaviour Atrey et al. [2019].

Our main contributions are:

  * •

We introduce ATOMs that quantify an agent’s attention to objects in its
environment, using _Hierarchical-attention_ and _Combinatorial-attention_
metrics.

  * •

We present three variation of the Pong game, each designed to induce distinct
behaviour in agents.

  * •

We show that ATOMs can provide insights into what an agent has learned about
its game and that successfully trained agents developed similar, game-
specific, attention patterns.

  * •

We used ATOMs to follow the development of an agent’s attention during its
training and identified distinct phases of attention patterns common to all
agents in our experiments. Additionally, we found that the agent’s self-
attention developed late in the learning and coincided with a significant
increase in its performance score.

##  2 Related work

#### Explainable RL

Explainable RL methods have been developed to elucidate the decision-making
process in RL agents. Some methods aim to enhance the interpretability of the
policy an agent has learned. This is achieved by encoding the policy in a
human-understandable form, such as a decision tree Liu et al. [2019], Coppens
et al. [2019], Roth et al. [2019], a set of linguistic rules Hein et al.
[2017] or a domain specific (programming) language Verma et al. [2018]. Other
methods decompose a complex policy into several simpler ones such as in Shu et
al. [2017] where the general policy is decomposed into local policies
corresponding to elementary skills that have been manually annotated by
humans. The global, more complex policy is then employed to determine which
local policy to apply along the agent’s trajectory. Alternative techniques
propose a justification for the model actions using structural causal models
Madumal et al. [2020] or a justification for the general policy with a
language understandable by humans Hayes and Shah [2017], van der Waa et al.
[2018]. Last but not least, a number of methods rely on visual explanations of
an agent’s behaviour. These include projection methods Zahavy et al. [2016]
that use the agent’s neural population activity to display the pattern learned
by its network. Whilst previous studies on saliency maps tend to focus on
fully trained models, we propose using them to study the learning dynamics
that occur during model training.

#### Saliency maps

Saliency maps have widely been used to explain predictions of deep image
classifiers in computer vision, but have been increasingly used as an
exploratory tool for RL models Zeiler and Fergus [2014], Zhou et al. [2016],
Weitkamp et al. [2018], Selvaraju et al. [2017]. There are two main types of
saliency maps: gradient based Greydanus et al. [2018] and perturbation-based
Puri et al. [2019] methods. Gradient-based methods estimate the significance
of each input pixel on the network’s decision-making process by analyzing the
gradient with respect to these pixels. In contrast, perturbation-based methods
estimate this influence by individually altering each pixel and observing the
resultant effect on the network’s choice of action. Nonetheless, perturbation-
based methods necessitate multiple forward passes for each input, which rend
them significantly computationally expensive Xing et al. [2023]. Gradient-
based methods are faster, but are also subject to noise. We chose Layer-wise
Relevance Propagation (LRP) Bach et al. [2015], Lapuschkin et al. [2019], as
it can be computed efficiently, while overcoming the noise problem many
gradient methods share.

##  3 Methods

###  3.1 Experimental setup

Data were recorded in actor-critic agents (A2C; Konda and Tsitsiklis [1999])
from the Stable Baseline 3 repository Raffin et al. [2021]. We study the actor
network of each agent, responsible for the choice of action. This network is
composed of three convolutional layers (size 4x32, 32x64 64x64) followed by a
linear layer (size 1x512) that we will refer to as $F_{c}$ (Fig 1 a) and a
final output layer. Throughout the remainder of this paper, we will use the
terms ’agent’ and ’actor network’ synonymously. The agent’s strategy was
derived from the activity of neurons in the linear layer ($F_{c}$). This layer
is located just before the output layer and encompasses the final world-model
on which the agent choose its action. ATOMs are computed from this layer.

###  3.2 Attention-Oriented Metrics (ATOMs)

ATOMs are composed of two metrics evaluating different aspects of an agent’s
attention: the hierarchical-attention and the combinatorial-attention metrics.
The hierarchical-attention evaluates the degree of attention given to an
individual object compared to the other objects. In contrast, the
combinatorial-attention evaluates how much attention an agent pays
simultaneously to different combinations of objects. This simultaneous
attention could be an indicator of a latent variable encoded by the neuron.
Here we characterize the Pong game by its constituent objects: the paddles of
the agent and its opponent, their respective displayed score, both balls (B1
and B2) and the walls. Figure 1 a illustrates the process to extract these two
metrics.

Both metrics composing ATOMs are derived from the attention of the neurons in
the $F_{c}$ layer that are relevant to an action. Specifically, this attention
is measured by computing Layer-Wise Relevance Propagation (LRP; Bach et al.
[2015], Lapuschkin et al. [2019]). For a model output $f(\boldsymbol{x})$
(here at the action node), LRP calculates a relevance score for each neuron in
the previous layer. Due to its conservation property, relevance can be
propagated layer-by-layer up until the input (here images) resulting in a
relevance map per layer. To generate input images $\boldsymbol{x}$ (consisting
of pixels p) used for the LRP analysis, we let an agent play the game for ten
episodes. Subsequently, we filtered the input to ensure all objects are
present and have a spacing of 10 pixels between them to prevent overlap. We
then sampled a total of 200 frames and automatically labeled the objects
within each frame. Let $\boldsymbol{X}\in\mathbb{R}^{200\times 4\times
84\times 84}$ and $\boldsymbol{\bar{X}}\in\mathbb{R}^{200\times 4\times
84\times 84}$ be the original and the labeled version of the input set,
respectively. We focus on the attention of neurons influencing action choices,
identifying relevant neurons in the $F_{c}$ layer using LRP.

For each input, $\boldsymbol{x}\in\boldsymbol{X}$:

  1. 1.

Calculate the relevance scores $R_{F_{c}}=\\{R_{i}\\}_{i\in F_{c}}$ and select
neurons $S\subseteq F_{c}$ that account for 90% of the total relevance
$R_{F_{c}}$. This filtering step is used to avoid noise components and
decrease ATOMs computational time.

  2. 2.

Generate relevance maps at the model input
$\boldsymbol{R}^{k}(\boldsymbol{x})=\\{R_{p}^{k}(\boldsymbol{x})\\}$ for each
neuron $k\in S$

  3. 3.

Identify the objects highlighted by the relevance maps.

#### Hierarchical-attention

We first characterize an agent’s attention by identifying which objects are
most important to the network’s decision-making. We do this by using a
hierarchical-attention metric $h:O\rightarrow\mathbb{R}$, to measure the
importance of each object. For objects $o_{g}\in O$, the hierarchical-
attention metric $h$ and its components are defined as follows:

$\displaystyle h(o_{g})$
$\displaystyle=\frac{1}{|\boldsymbol{X}|}\sum_{\boldsymbol{x}\in\boldsymbol{X}}\sum_{k\in
S}R_{k}(\boldsymbol{x})\cdot\bar{R}_{\text{g}}^{k}(\boldsymbol{x}),$ (1)
$\displaystyle\bar{R}_{\text{g}}^{k}(\boldsymbol{x})$
$\displaystyle=\frac{1}{V}\sum_{p\in
P_{o_{g}}}R_{p}^{k}(\boldsymbol{x}),\quad\text{where
}V=|R_{p_{g}}^{k}(\boldsymbol{x})\setminus\\{0\\}|,$ (2)

and $p_{g}$ denotes all pixels $p\in o_{g}$.

#### Combinatorial-attention

We then incorporate the combinations of objects that an agent simultaneously
attends to by using the combinatorial-attention metric
$c:O\rightarrow\mathbb{R}$. For a subset of objects $T\subseteq O$, the
combinatorial-attention metric is defined as follows:

$\displaystyle c(T)$
$\displaystyle=\frac{1}{|\boldsymbol{X}|}\sum_{\boldsymbol{x}\in\boldsymbol{X}}\sum_{k\in
S}R_{k}(\boldsymbol{x})\cdot\delta_{k}(\boldsymbol{x};T),$ (3)
$\displaystyle\delta_{k}(\boldsymbol{x};T)$
$\displaystyle=\begin{cases}1&\text{if
}\bar{R}_{\text{g}}^{k}(\boldsymbol{x})>\beta\text{ for all }g\in T\text{ and
}\bar{R}_{\text{g}}^{k}(\boldsymbol{x})=0\text{ for all }g\in O\setminus T,\\\
0&\text{else.}\end{cases}$ (4)

where $\beta=\alpha\cdot M^{k}$ with $\alpha=0.25$ and
$M^{k}=\max_{g}(\\{\bar{R}_{\text{g}}^{k}(\boldsymbol{x})\\}_{g})$.

An illustration of different combination of objects looked at by neurons in
the $F_{c}$ layer can be found in appendix
LABEL:combinatorial_attention_example (figure
LABEL:fig:example_combin_attention). The impact of the $\alpha$ threshold on
the results of the combinatorial-attention metric is discussed in appendix
LABEL:threhold_value_combinatorial_attention.

###  3.3 Measurements of ATOMs during training

The frequency of evaluation ($f_{r}$), denoting how often the metrics were
extracted, was chosen for each game based on the training dynamics observed.
We set $f_{r}=230$ for version $v_{0}$, $f_{r}=85$ for version $v_{1}$ and
$f_{r}=50$ for version $v_{2}$. The choice of $f_{r}$ for each game was
influenced by the number of epochs required for the network to demonstrate an
increase in training score in the respective game. It’s important to note that
different frequencies could be selected based on the desired granularity for
studying the evolution of the metrics. Our selected frequencies aimed to
balance the need for detailed metric tracking with the practical constraints
of computational resources and training time.

###  3.4 Designing Pong variations

We implemented three different variations of the Pong game: v0, v1, and v2.
Each variant builds on the classic Pong gameplay where two agents compete with
each other, each controlling a paddle to score points by returning the ball to
the opponent’s side. A point is awarded when the ball passes the opponent’s
paddle and an episode terminate when a player accumulates 21 points. For
training, a reward of +1 or -1 is assigned depending on which side scores. v0
maintains traditional gameplay. In v1, two balls are present (B1 and B2), but
only B1 counts for scoring while B2 is used as a distraction. In v2, both B1
and B2 are integral, but B2 only interacts with the agent’s paddle, awarding
+1 for a rebound and -1 if it passes the paddle. Additionally, in v2, the
agent must reach 41 points, while the opponent needs only 21 to win. Figure 1b
illustrates the three game variations. Details regarding the game
implementation and the rules can be found in Appendix LABEL:Game.

Figure 1: a ATOMs consists of two metrics (hierarchical, and combinatorial
attention) that summarize the agent’s attention. These metrics are derived
from relevance maps generated using the LRP method, applied to a typical game
state, where all objects are simultaneously presented. b Illustration of the
three variations of the Pong game. Red balls yield reward to the agent. In v1
the gray ball is a distractor. The balls’ dynamics are illustrated with lines
indicating if a ball bounces back from the opponent or passes through it.

###  3.5 Evaluating agent behaviour: Dual Ball Discrimination Test

A major difference between the three game variations is the importance of ball
B1 and ball B2. As the latter only yield rewards in v2 we expect the behaviour
of the agents toward B2 to change between games. To quantitatively evaluate
agent behaviour, we implemented the Dual Ball Discrimination Test, wherein the
agent is put in a situation enforcing a choice between both balls. This
experiment involved generating 100 unique trajectories from various initial
positions and velocities, designed such that both balls reach the agent’s
x-coordinate simultaneously while maintaining a separation exceeding the
length of the paddle. Consequently, the agent can interact with only one of
the balls and we can analyse any systematic preference for B1 or B2. An
illustration of the test is displayed figure 3 a.

##  4 Experiments

###  4.1 Evaluation of ATOMs

Figure 2: a Hierarchical profile averaged over 50 trained agents for each game
variation. b Dissimilarity matrix for the hierarchy-attention metrics for all
agents.c Combinatorial-attention metric which examines co-observation of
objects by individual neurons. Each object is symbolized by a distinct
coloured dot. Combinations of objects are indicated through the simultaneous
colouring of two or more dots. d Dissimilarity matrix for the combinatorial-
attention metrics for all agents. e Dissimilarity between agents within game
variation for the hierarchy and the combinatorial-attention metrics. Error
bars represent standard-deviation values.

We evaluated the capacity of ATOMs to provide relevant information about what
an agent has learned. To this end, we computed ATOMs on 50 fully trained
agents for each game variation: v0, v1, v2. Agents were trained using
different random initialization (details regarding agent training protocols
are provided in LABEL:RL-training-procedure).

#### Different game variations induce game-specific attention pattern.

The hierarchical-attention (figure 2 a) shows that trained agents attended to
the expected aspects of the game. Notably, agents trained on v0 and v1 payed
more attention to ball B1 (the only ball that brings a reward) while agents
trained on v2 (in which both balls yield rewards) showed a more nuanced
attention pattern with a slight preference for the ball B2. Across all game
variations, the second highest attention was consistently directed towards the
agent and opponent paddles. Interestingly, agents from v2 developed attention
to the score displayed on the game screen as well.

The combinatorial-attention (figure 2 c) further breaks down the attention of
the neurons in the $F_{c}$ layer into combination of objects. We found a
consistent common pattern across game variations with attention to either the
ball(s) alone or to a combination of ball(s) and the agent’s paddle or ball(s)
and the opponent’s paddle. Figures 2 b,d show that successfully trained agent
developed similar attention patterns as the similarity in both ATOMs metrics
are higher within game variation. Interestingly, we observed greater
variability in the combinatorial-attention metrics (figure 2 e) in variation
v2 compared to v1 and v0 indicating that agents may develop more diverse
attention patterns in game variations that are more complex.

#### Attention patterns are consistent with agents’ behaviour.

As shown in figure 2 a, ATOMs indicate that agents trained on v1 still paid
attention to B2 suggesting that these agents had not learned to ignore the
ball completely. To investigate this further, we examined whether the
attention directed towards B1 and B2 reflected the agent’s behaviour, which we
defined here as the relative interaction of the agent with each of the balls.
Specifically, we assessed whether increased attention to one ball over the
other reflected a preference for interacting with that particular ball. To
study this experimentally, we exposed agents trained on game variation v1 or
v2 to the Dual Ball Discrimination Test 3.5 to measure how often the agent
interacted with each ball (figure 3 a). Figure 3 b shows he proportion of B1
interactions relative to B2 interactions against the importance given to B1
relative to B2. Distinct clusters emerge that correspond to variations v1 and
v2. This shows that the relative attention given to each ball indeed reflects
how an agent behaves towards them (in terms of preferred interaction).
Furthermore, the cluster associated with v1 shows that these agents still
interact with B2, confirming the hypothesis driven by the ATOMs metrics that
agents in v1 had not learned to ignore B1 completely.

As mentioned, agents trained on v2 developed an attention to the score
displayed on the screen (figure 2 a). This prompted us to analyse the impact
of the displayed score on the agent’s behaviour. While this may be an artifact
of an increased spatial attention in the more complex game variant, the
results of the  Dual Ball Discrimination Test reveal that the displayed scores
indeed influenced the agent’s choice to hit B1 or B2. Distinct patterns of
influence were observed among agents. For instance, some agents exhibited a
preference for B1 until their displayed score approached 40 points, at which
point they tended to choose B2 (the game terminates when the agent reaches 41
points). Other agents shifted their preference from B1 to B2 once the
displayed score of the opponent exceeded 10 points. Further details of the
experiments and results are provided in Appendix LABEL:impact_score_choice_v2.

Figure 3: a Dual Ball Discrimination Test, 100 trajectories were generated
such that the agent is forced to choose between interacting with B1 or B2. b
Relative interaction of B1 with respect to B2 in function of the relative
hierarchy of B1 with respect to B2. v0 was added as a reference.

###  4.2 Evolution of agents’ attention during learning

We then monitored the development of the agents’ attention during their
training using ATOMs.

Figure 4:  Development of the performance score and ATOMs during learning for
each game version. The performance scores of the agent were computed for ten
games at each point of measurement during training. The average score is
depicted by the central dark line, with the shaded region representing $+/-$
the standard deviation. The dotted lines indicate the time frame during which
the agent’s score showed a marked increase. These time frames were picked
manually. The combinatorial-attention is categorized as follows from top to
bottom: noise (grey), combinations of balls only (orange), combinations of
opponent and balls (cyan),and combinations of objects including the agent
itself (red)

.

#### Application of ATOMs during training

Figure 4 illustrates ATOMs for three agents trained on v0, v1 and v2 alongside
their performance during training. The results reveals that most pronounced
changes in the hierarchical metric occurred at the initial stage of the
training preceding a marked enhancement in the agent’s performance. While the
hierarchical metric indicates the overall importance of an object, the
combinatorial metrics dissect this information further by revealing to which
combination of objects an agent pays attention to. In this figure, we
categorized the plots of the combinatorial metrics into four major groups:

  * •

combinations involving only B1 and B2 (Ball(s)),

  * •

combinations involving the balls and the opponent (Opponent),

  * •

combinations involving the agent (Agent),

  * •

combinations of objects only present during the initial learning phase
(Noise).

The evolution of the combinatorial-attention during training shows how the
attention of an agent gets progressively selective to specific subsets of
objects. For example, by looking at the evolution of the attention to the
ball(s) for the games v1 and v2, one can see that the agent’s attention
started by including both ball objects (B1, B2). Then as the training progress
its attention progressively isolated the individual balls (B1 in v1 and B1 and
B2 in v2).

#### Attention patterns develop in phases during learning

The four categories defined above appeared to follow a timeline consistent
amongst versions v0, v1 and v2 (Figure 4 ). We will therefore refer to these
categories as phases. In the first phase the attention was non-selective and
spread among objects (Noise phase). In the following phase the agents focused
most attention upon the moving ball(s). Subsequently, the agents also paid
attention to the opponent. Finally, in the "Agent phase", agents payed
attention to themselves in relation to the other objects. Interestingly, this
phase was consistently the last one to emerge and coincide with a marked
increase in the agents’ performance scores (indicated by the two dotted lines
in figure 5) 111Agents trained on game v2 developed an additional attention
component that integrated the scores of the agent and/or of the opponent
(results in appendix LABEL:fig:learning_v2). To test if theses phases were
consistent among all agents, we repeated this analysis with ten agents per
version and measured their phase onsets ( figure 5). We defined the onset of a
phase as the point at which the first component composing this phase exceeds
25% of its maximum value. This onset was then expressed as a percentage of the
total training period for a given agent, with the total training period being
the time taken for an agent to reach its maximum score. We found that the
timeline of these attention phases were consistent across models and game
versions.

Figure 5:  Onset of the attention phases (Noise, Ball(s), Opponent, Agent) for
each agent (represented by different colours) in game versions textit{v0},
textit{v1}, and textit{v2}. Each phase is composed of one or several
components . The onset of a phase is defined as the point at which the first
component composing this phase exceeds 25% of its maximum value. This onset is
then expressed as a percentage of the total training period for a given agent.

##  5 Discussion

The learning process of an RL agent is poorly understood beyond the
mathematical equations of its learning algorithm. In this study we delve
deeper into the learning phase of an agent by analyzing the development of its
attention. To this end, we introduce Attention-Oriented Metrics (ATOMs)
designed to quantitatively assess an agent’s attention. We evaluated ATOMs’
ability to give information about what an RL agent has learned, by
implementing three versions of a Pong game, each requiring the agent to learn
distinct behaviour.

Our results demonstrated that ATOMs successfully delineate the different
attention patterns that each game version induces. More specifically, we
observed that successfully trained agents exhibited consistent, game-specific
attention patterns. Additionally, behavioural tests revealed that these
divergent attention patterns corresponded to observable differences in agent
behaviour. Through continuous monitoring of ATOMs during training, we
identified distinct, sequentially ordered developmental phases in its
attention. Notably, an agent’s attention to itself emerged relatively late in
the learning process coinciding with an increase in its performance score. We
suggest that an agent’s attention is guided by the association of objects’
movements - or more generally changes in the state of an object- with reward
feedback. Indeed, in our three games, the most straightforward correlation
between an object’s movement and a feedback reward involves the movement of
the ball(s), either by passing the agent side or the opponent’s side, or by
hitting the agent’s paddle (for B2 in version v2). Since the opponent is
designed to score points, its movement is correlated with future rewards as
well. Conversely, as the agent learns to play pong its actions are less likely
to influence the feedback reward. This difference could explain the agent’s
initial focus on the ball(s) and the opponent and the late focus to the paddle
it controls.

While ATOMs offer valuable insights about RL agent development during
learning, they have limitations that could be addressed in future work.
Indeed, ATOMs are object-based, which, while reflecting our human
understanding of the games, may limit their explanatory power. Alternative
metrics based on an agent’s attention, not limited to the object level, could
reveal more subtle abstractions or reasoning underlying an agents’ decision-
making processes. Additionally, ATOMs requires image labeling and are best
suited to environments with a limited set of objects, as its combinatorial-
attention component grows polynomially with the number of objects potentially
diminishing its explanatory power.

Despite these constraints, we believe that ATOMs can enhance our understanding
of the learning process of RL agents, and offer several potential
applications. For instance, monitoring an agent’s attention could ensure its
focus on the pertinent objects (either during training or online) particularly
in safety-critical applications. It could also facilitate the early detection
of flawed policies or exploitation of bugs in an environment which might not
be evident through performance scores alone. For example, in our experiments,
agents trained on v1 had not learn to fully ignore the irrelevant ball despite
achieving an optimal performance score. This suggest potential refinement in
the environment design or in the reward function. Thus, accessing an agent’
attention during its training could accelerate the debugging of the
environments and reward functions.

## Acknowledgments

N.S., C.B. and S.H. are supported by BMBF (Federal Ministry of Education and
Research) through ACONITE (01IS22065) and the Center for Scalable Data
Analytics and Artificial Intelligence (ScaDS.AI.) Dresden/Leipzig. C.B. is
also supported by the Max Planck IMPRS NeuroCom Doctoral Program. The figures
were created using BioRender.com.

## References

  * Watkins and Dayan [1992] Christopher JCH Watkins and Peter Dayan.  Q-learning.  _Machine learning_ , 8:279–292, 1992.
  * Konda and Tsitsiklis [1999] Vijay Konda and John Tsitsiklis.  Actor-critic algorithms.  _Advances in neural information processing systems_ , 12, 1999.
  * Mnih et al. [2013] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Alex Graves, Ioannis Antonoglou, Daan Wierstra, and Martin Riedmiller.  Playing atari with deep reinforcement learning.  _arXiv preprint arXiv:1312.5602_ , 2013.
  * Clark and Amodei [2016] Jack Clark and Dario Amodei.  Faulty reward functions in the wild.  _Internet: https://blog. openai. com/faulty-reward-functions_ , 2016.
  * Cobbe et al. [2019] Karl Cobbe, Oleg Klimov, Chris Hesse, Taehoon Kim, and John Schulman.  Quantifying generalization in reinforcement learning.  In _International Conference on Machine Learning_ , pages 1282–1289. PMLR, 2019.
  * Verma et al. [2018] Abhinav Verma, Vijayaraghavan Murali, Rishabh Singh, Pushmeet Kohli, and Swarat Chaudhuri.  Programmatically interpretable reinforcement learning.  In _International Conference on Machine Learning_ , pages 5045–5054. PMLR, 2018.
  * Shu et al. [2017] Tianmin Shu, Caiming Xiong, and Richard Socher.  Hierarchical and interpretable skill acquisition in multi-task reinforcement learning.  _arXiv preprint arXiv:1712.07294_ , 2017.
  * Liu et al. [2019] Guiliang Liu, Oliver Schulte, Wang Zhu, and Qingcan Li.  Toward interpretable deep reinforcement learning with linear model u-trees.  In _Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2018, Dublin, Ireland, September 10–14, 2018, Proceedings, Part II 18_ , pages 414–429. Springer, 2019.
  * Madumal et al. [2020] Prashan Madumal, Tim Miller, Liz Sonenberg, and Frank Vetere.  Explainable reinforcement learning through a causal lens.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 34, pages 2493–2500, 2020.
  * Hein et al. [2017] Daniel Hein, Alexander Hentschel, Thomas Runkler, and Steffen Udluft.  Particle swarm optimization for generating interpretable fuzzy reinforcement learning policies.  _Engineering Applications of Artificial Intelligence_ , 65:87–98, 2017.
  * Zeiler and Fergus [2014] Matthew D Zeiler and Rob Fergus.  Visualizing and understanding convolutional networks.  In _European conference on computer vision_ , pages 818–833. Springer, 2014.
  * Zhou et al. [2016] Bolei Zhou, Aditya Khosla, Agata Lapedriza, Aude Oliva, and Antonio Torralba.  Learning deep features for discriminative localization.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 2921–2929, 2016.
  * Weitkamp et al. [2018] Laurens Weitkamp, Elise van der Pol, and Zeynep Akata.  Visual rationalizations in deep reinforcement learning for atari games.  In _Benelux Conference on Artificial Intelligence_ , pages 151–165. Springer, 2018.
  * Selvaraju et al. [2017] Ramprasaath R Selvaraju, Michael Cogswell, Abhishek Das, Ramakrishna Vedantam, Devi Parikh, and Dhruv Batra.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  In _Proceedings of the IEEE international conference on computer vision_ , pages 618–626, 2017.
  * Atrey et al. [2019] Akanksha Atrey, Kaleigh Clary, and David Jensen.  Exploratory not explanatory: Counterfactual analysis of saliency maps for deep reinforcement learning.  _arXiv preprint arXiv:1912.05743_ , 2019.
  * Coppens et al. [2019] Youri Coppens, Kyriakos Efthymiadis, Tom Lenaerts, Ann Nowé, Tim Miller, Rosina Weber, and Daniele Magazzeni.  Distilling deep reinforcement learning policies in soft decision trees.  In _Proceedings of the IJCAI 2019 workshop on explainable artificial intelligence_ , pages 1–6, 2019.
  * Roth et al. [2019] Aaron M Roth, Nicholay Topin, Pooyan Jamshidi, and Manuela Veloso.  Conservative q-improvement: Reinforcement learning for an interpretable decision-tree policy.  _arXiv preprint arXiv:1907.01180_ , 2019.
  * Hayes and Shah [2017] Bradley Hayes and Julie A Shah.  Improving robot controller transparency through autonomous policy explanation.  In _Proceedings of the 2017 ACM/IEEE international conference on human-robot interaction_ , pages 303–312, 2017.
  * van der Waa et al. [2018] Jasper van der Waa, Jurriaan van Diggelen, Karel van den Bosch, and Mark Neerincx.  Contrastive explanations for reinforcement learning in terms of expected consequences.  _arXiv preprint arXiv:1807.08706_ , 2018.
  * Zahavy et al. [2016] Tom Zahavy, Nir Ben-Zrihem, and Shie Mannor.  Graying the black box: Understanding dqns.  In _International conference on machine learning_ , pages 1899–1908. PMLR, 2016.
  * Greydanus et al. [2018] Samuel Greydanus, Anurag Koul, Jonathan Dodge, and Alan Fern.  Visualizing and understanding atari agents.  In _International conference on machine learning_ , pages 1792–1801. PMLR, 2018.
  * Puri et al. [2019] Nikaash Puri, Sukriti Verma, Piyush Gupta, Dhruv Kayastha, Shripad Deshmukh, Balaji Krishnamurthy, and Sameer Singh.  Explain your move: Understanding agent actions using specific and relevant feature attribution.  _arXiv preprint arXiv:1912.12191_ , 2019.
  * Xing et al. [2023] Jinwei Xing, Takashi Nagata, Xinyun Zou, Emre Neftci, and Jeffrey L Krichmar.  Achieving efficient interpretability of reinforcement learning via policy distillation and selective input gradient regularization.  _Neural Networks_ , 161:228–241, 2023.
  * Bach et al. [2015] Sebastian Bach, Alexander Binder, Grégoire Montavon, Frederick Klauschen, Klaus-Robert Müller, and Wojciech Samek.  On pixel-wise explanations for non-linear classifier decisions by layer-wise relevance propagation.  _PloS one_ , 10(7):e0130140, 2015.
  * Lapuschkin et al. [2019] Sebastian Lapuschkin, Stephan Wäldchen, Alexander Binder, Grégoire Montavon, Wojciech Samek, and Klaus-Robert Müller.  Unmasking clever hans predictors and assessing what machines really learn.  _Nature communications_ , 10(1):1096, 2019.
  * Raffin et al. [2021] Antonin Raffin, Ashley Hill, Adam Gleave, Anssi Kanervisto, Maximilian Ernestus, and Noah Dormann.  Stable-baselines3: Reliable reinforcement learning implementations.  _Journal of Machine Learning Research_ , 22(268):1–8, 2021.  URL http://jmlr.org/papers/v22/20-1364.html.