# Resistive Hose Modes in Tokamak Runaway Electron Beams

A.P. Sainterme<EMAIL_ADDRESS>C.R. Sovinec  Department of Nuclear
Engineering & Engineering Physics, University of Wisconsin - Madison, 1500
Engineering Drive, Madison, Wisconsin 53706-1609

###### Abstract

Beams of energetic runaway electrons are generated during disruptions in
tokamaks, and fluid models are used to study their effects on macroscale
dynamics. Linear computations of a massless, runaway electron beam coupled to
MHD plasma show that resistive hose instabilities grow faster than tearing
modes at large resistivity. Eigenvalue results with reduced models of the
resistive hose instability are compared with results from the full MHD and
beam system, showing that the resistive hose decouples from any plasma
response. An estimate of plasma temperature at which growth of the resistive
hose dominates tearing for post-disruption DIII-D plasma parameters is in a
physically relevant regime

The tokamak is the leading candidate for a magnetic confinement fusion
reactor. It confines plasma in a toroidal configuration with an externally
generated toroidal magnetic field and a poloidal magnetic field supported by
toroidal electrical current that is driven through the plasma. Tokamaks are
prone to disruptions—sudden dynamic events in which confinement of the plasma
is lost and the plasma current is terminated. A typical disruption begins with
a rapid loss of plasma thermal energy, a thermal quench (TQ), followed by a
current quench (CQ). The decrease in temperature during the TQ increases the
electrical resistivity of the plasma, increasing the electric field along
magnetic field-lines and accelerating a small population of electrons to high
energies. Since the effective collision frequency decreases at higher
energies, the resulting beam-like population of supra-thermal electrons is
largely unimpeded by collisional drag. These runaway electrons (REs) can carry
a substantial fraction of the initial plasma current and have caused
considerable damage to plasma-facing componentsBreizman  _et al._ (2019).
Understanding the dynamics of REs during disruptions is critical for the
development of tokamak fusion reactors. For the ITER experiment, which will
carry a plasma current of 15 MA, being able mitigate RE-induced damage will be
essentialLehnen  _et al._ (2015). To that end, the development and
characterization of computational tools to simulate RE dynamics in tokamaks is
required.

One common method for modeling the macroscale effects of REs in tokamaks
treats them as a separate cold beam-like fluid species in extended MHD
simulations. The RE beam is treated as a source of resistance-free current
density whose direction depends on the time-evolving magnetic field that
interacts with a plasma consisting of thermal electrons and ions governed by
single-fluid MHD equations. The consistency of this approximation presumes
negligible inertia of both electron species, $m_{e}\ll m_{i}$,
$\gamma_{r}m_{e}\ll m_{i}$, and a small particle density of REs
$n_{r}\ll\rho/m_{i}$. Here, $\gamma_{r}$ is the relativistic factor associated
with the runaway parallel velocity: $\gamma_{r}=1/\sqrt{1-(c_{r}/c)^{2}}$.
Since the RE current is essentially resistance freeBreizman  _et al._ (2019),
the typical resistive MHD Ohm’s law is modified:

$\mathbf{E}+\mathbf{V}\times\mathbf{B}=\eta\mathbf{J}\;\rightarrow\;\mathbf{E}+\mathbf{V}\times\mathbf{B}=\eta\left(\mathbf{J}-\mathbf{J}_{r}\right).$

The runaway current density is subtracted from the total current density so
that only the fraction of current carried by the bulk plasma is affected by
resistivity.

Because the resistivity is modified, a natural question is the effect of the
inclusion of REs on resistive MHD instabilities. Prior work with this model
has explored the effect of a fluid runaway current on the linear behavior of
the tearing mode in a sheet pinch Avinash and Kaw (1988), and tearing and
resistive kink modes in a large aspect ratio cylinder Liu  _et al._ (2020)
Zhao  _et al._ (2020). Helander et al. additionally consider nonlinear
saturation of the tearing modeHelander  _et al._ (2007).

The stability of high-energy, relativistic particle beams propagating through
resistive, neutralizing background plasma has also been studied. The earliest
analytical description of an instability in such a system is due to Rosenbluth
Rosenbluth (1960). Rosenbluth identifies a kink-type instability of a self-
pinched beam of relativistic particles that couples the motion of the beam
particles with the resistive response of magnetic field. Weinberg generalizes
the work of Rosenbluth first to a spatially modulated beam Weinberg (1964),
where the instability is dubbed a ‘hose’ mode, and later, to a more general
dispersion relation for the self-pinched equilibrium Weinberg (1967) The hose
mode has long been recognized in the accelerator community as an impediment to
the propagation of beams, and controlling its growth with properly tailored
beam profiles has been a subject of studyFernster  _et al._ (1995).

In the resistive hose case, it is assumed that the plasma column surrounding
the beam is generated by ionization of a dense gas by the beam particles. The
beam is self-pinched, meaning that the self-induced magnetic field is
azimuthal, and the equilibrium force balance condition requires that this
field is sufficiently strong to provide the centripetal acceleration of the
beam particle orbits around the axis. Some relatively recent work includes the
effect of an axial magnetic field, which modifies the equilibrium orbits of
the beam particlesUhm and Lampe (1980). In contrast, in the case of the RE
beam generated in a tokamak, the background plasma carries current density
before the RE population is generated. There is an existing MHD equilibrium
with both axial (toroidal) and azimuthal (poloidal) magnetic field, and after
its generation, the RE beam carries some substantial fraction of the
equilibrium current density. Despite their differences, we show that important
aspects of the resulting mathematical descriptions of the particle-beam and
tokamak RE systems are equivalent.

In this letter, we draw the connection between the earlier results on
resistive hose instabilities of particle beams and resistive instabilities in
the context of RE modeling in tokamaks. It is notable that the dynamics of the
RE beam in this work are treated under the assumption that the electrons have
negligible mass and are guided by the tokamak magnetic fields in contrast to
the particle-beam analysis, where inertia is important for the particle
orbits. No analysis of the resistive hose in the literature has, heretofore,
explored the massless limit. To check its relevance, we estimate that the
resistive hose mode that appears in this work would grow faster than the
tearing mode at post-TQ plasma temperatures of $T_{e}\approx 1.5$ eV or less
in DIII-D Lvovskiy  _et al._ (2023).

The model consists of a background plasma whose dynamics are governed by
single-fluid, resistive MHD. In the present discussion, the effects of
pressure are presumed to be negligible, as appropriate for post-TQ conditions.
The MHD equations are augmented with a second fluid species that represents
the collective low-frequency behavior of REs. The model for REs is an electron
beam with a large constant speed parallel to the magnetic field. In the
present work, the $\mathbf{E}\times\mathbf{B}$ and curvature drifts are
neglected. The RE beam dynamics are coupled to the resistive MHD equations via
the modification to the Ohm’s law noted earlier. Quasi-neutrality between the
species is assumed so that there are no space-charge effects from the beam.

The dynamical equations are linearized about a steady-state solution, and
equations are derived for the time-evolution of small perturbations from the
equilibrium.

$\displaystyle\partial_{t}n_{r}+\nabla\cdot\left(n_{r}\mathbf{U}_{r}+N_{r}\mathbf{u}_{r}\right)=0,$
(1)
$\displaystyle\rho\partial_{t}\mathbf{v}=\mathbf{j}\times\mathbf{B}+\mathbf{J}\times\mathbf{b},$
(2) $\displaystyle\partial_{t}\mathbf{b}=-\nabla\times\mathbf{e},$ (3)
$\displaystyle\nabla\times\mathbf{b}=\mathbf{j},$ (4)
$\displaystyle\mathbf{e}=-\mathbf{v}\times\mathbf{B}+\eta\left(\mathbf{j}-\mathbf{j}_{r}\right),$
(5) $\displaystyle\mathbf{U}_{r}=-c_{r}\mathbf{B}/B,$ (6)
$\displaystyle\mathbf{u}_{r}=-c_{r}\mathbf{b}_{\perp}/B,$ (7)
$\displaystyle\mathbf{j}_{r}=-e\left(n_{r}\mathbf{U}_{r}+N_{r}\mathbf{u_{r}}\right),$
(8) $\displaystyle\nabla\cdot\mathbf{b}=0.$ (9)

Equations 1-9 describe the time dependence of the perturbations to the
magnetic field, $\mathbf{b}$, the bulk plasma velocity, $\mathbf{v}$, and the
runaway density, $n_{r}$. The variables with capital letters and the mass
density $\rho$ represent the equilibrium fields, and the lowercase variables
are the perturbations.

As mentioned previously, the impetus for this model is to capture the effect
that a presence of a substantial amount of runaway electron current has on
instabilities that are present in resistive MHD. An early analytic application
of this model to tearing modes in slab geometry found that in the limit of
zero inertia, stability of the tearing mode is unaffectedAvinash and Kaw
(1988). More recent results using reduced MHD in cylindrical geometry also
suggest that the runaway current does not impact the stability of the tearing
mode, but that the nonlinear saturation is affectedHelander  _et al._ (2007).
The analytic and numerical work in Liu  _et al._ (2020) finds an additional
correction to the linear dispersion relation that results in overstability of
the tearing mode.

We have implemented equations 1-9 in the NIMROD extended MHD codeSovinec  _et
al._ (2004). The linear equations are Fourier analyzed in the axial direction,
and a high-order spectral element representation is used in the poloidal
$(r-\theta)$ plane. For a single axial Fourier harmonic, $e^{ikz}$, the
initial value problem is solved by numerically integrating the linear
equations in time from an arbitrary initial condition. As time
$t\rightarrow\infty$, the solution approaches the most unstable mode from
which the growth rate and frequency are determined.

For cylindrical problems, the equations can be reduced to a set of coupled
ordinary differential equations in the radial coordinate by also assuming a
Fourier representation in the azimuthal angle, $e^{im\theta}$ and a complex
exponential dependence in time $e^{\gamma t}e^{-i\omega t}$. To verify NIMROD
results, a 1D eigenvalue code is used to numerically solve the resulting
ordinary differential equations in $r$ to determine $\omega$ and $\gamma$
given $k$ and $m$. The eigenvalue solver also uses spectral elements, albeit
for the single spatial dimension $r$Sovinec (2016).

The plasma equilibrium considered here is a $0-\beta$ cylindrical screw pinch
with a peaked current density profile. The equilibrium current density is
assumed to be entirely carried by the RE beam. This equilibrium approximates
the situation observed tokamaks during the current plateau phase that follows
a TQHollmann  _et al._ (2013). The safety factor $q$, is given by
$q(r)=1.15(1+1.54(r/a)^{2})$. Note that $q(0)>1$, and that the profile is
monotonically increasing. The only MHD instabilities expected in the absence
of runaway effects are tearing modes with $m\geq 2$.

The cylinder is considered periodic in the axial coordinate $z$. It has unit
minor radius, $a=1$, there is a perfectly conducting wall placed at $r=a$, and
the RE beam radius is equal to the plasma radius. The axial field strength at
$r=0$ is $B_{0}=1$. The density $\rho$ is chosen so that the Alfvén velocity
$V_{A}=1$ when $B_{0}=1$, and $c_{r}=20V_{A}$. This particular value of
$c_{r}$ is chosen as a representative example for $c_{r}\gg V_{A}$. The
resistivity is normalized to $aV_{A}$. In the results that follow, the axial
wavenumber is $k=-0.1$.

Figure 1 plots the growth rate and frequency of the fastest growing mode from
NIMROD calculations of equations 1-9 as a function of varying resistivity
values. Also plotted are the growth rates of the fastest growing mode in the
absence of the RE beam, growth rates and frequencies computed with the M3D-C1
codeLiu  _et al._ (2020), and growth rates and frequencies computed with our
1D eigenvalue code with $m=2$.

Figure 1: Comparison of the growth rates and frequencies of the fastest
growing mode computed by NIMROD, M3D-C1, and the 1D eigenvalue code ($m=2$).
In the NIMROD result, the resistive hose mode is the fastest growing at
$\eta=10^{-3}$.

For values of $\eta\lesssim 10^{-4}$, the fastest growing mode is the $m=2$
tearing mode, which displays the expected localization around the $q=2$
surface. For $\eta>10^{-4}$, the fastest growing mode computed in NIMROD is a
hose-like mode with $m=1$. The appearance of this hose mode accounts for the
distinct growth rate and frequency at $\eta=10^{-3}$ reported by NIMROD.
Figure 2 plots contours of the radial and azimuthal components of the solution
for $\eta=10^{-2}$.

Figure 2: Contours of the radial and azimuthal components of the perturbed
magnetic field associated with the growing mode observed at $\eta=10^{-2}$.

The transition from the tearing mode to the hose mode is explicated by
calculating each $m$ separately with the eigenvalue code. Figure 3 plots the
frequencies and growth rates for $m=1$ and $m=2$ as a function of resistivity.
The shaded region denotes the range of $\eta$ values where the hose mode grows
faster than the tearing mode.

Figure 3: Growth rate and frequency of the fastest growing mode vs.
resistivity from the full MHD eigenvalue calculations for the $m=1$ beam mode
and the $m=2$ tearing mode.

It is clear from the eigenvalue code results that the hose mode scales
linearly with resistivity for small $\eta$, whereas the tearing mode scales as
$\eta^{3/5}$. The difference in scaling suggests that for equilibria that are
unstable to both modes, there is a value of $\eta$ at which the fastest
growing mode transitions from tearing to the hose. The exact details of the
transition depend on the equilibrium current profile, and quantifying this
will require an analytic expression for the hose dispersion relation in this
model. Deriving such an expression remains a task for future work.

Since the resistive hose mode is driven by the interaction of the beam with
the magnetic field, we can isolate the instability from the full system by
simply neglecting the MHD velocity perturbation $\mathbf{v}$. This is a good
approximation when the background plasma response is on a much longer time
scale than the hose mode dynamics. In making this simplification we also
redefine the RE density variables as

$\lambda\equiv\frac{ec_{r}n_{r}}{B},$ (10)

and

$\Lambda\equiv\frac{ec_{r}N_{r}}{B}.$ (11)

Then, when $\mathbf{v}=0$, we can consider only the equations

$\displaystyle{\partial_{t}\lambda-
c_{r}\frac{\mathbf{B}}{B}\cdot\nabla\lambda=\frac{c_{r}}{B}\nabla\cdot\left(\Lambda\mathbf{b}_{\perp}\right)},$
(12)
$\displaystyle{\partial_{t}\mathbf{b}-\eta\nabla^{2}\mathbf{b}=\eta\nabla\times\left(\lambda\mathbf{B}+\Lambda\mathbf{b}_{\perp}\right)}.$
(13)

Eigenvalues for equations 12 and 13 are computed with the same 1D eigenvalue
code. Note that in this system of equations, $V_{A}$ is not a relevant
parameter.

For the particular equilibrium considered, both with the full MHD system, and
with $\mathbf{v}=0$, the polarization of the perturbed magnetic field for the
hose-like mode is primarily transverse to the axis of the cylinder. i.e., the
axial component of the solution for the perturbed magnetic field is
negligible: $\mathbf{b}\cdot\hat{\mathbf{z}}\approx 0$. This observation
motivates the introduction of a Clebsch representation for the perturbed
magnetic field, $\mathbf{b}=\nabla\psi\times\hat{\mathbf{z}}$ . In this
approximation, again assuming the dependence $\psi=\psi(r)\exp(im\theta+ikz)$,
equations 12 and 13 reduce to two scalar equations for the complex quantities
$\psi$ and $\lambda$:

$\displaystyle\partial_{t}\psi-\eta\nabla^{2}\psi$ $\displaystyle=$
$\displaystyle\eta\left(B_{z}\lambda+\Lambda\frac{B_{\theta}B_{z}}{B^{2}}\psi^{\prime}\right)$
(14)
$\displaystyle+\eta\frac{kr}{m}\left(1-\frac{B_{\theta}^{2}}{B^{2}}\right)\psi^{\prime},$
$\displaystyle\partial_{t}\lambda-ic_{r}F\lambda$ $\displaystyle=$
$\displaystyle\frac{c_{r}}{B}\Lambda^{\prime}\frac{im}{r}\psi+ic_{r}F\Lambda\frac{B_{\theta}}{B^{2}}\psi^{\prime},$
(15)

where $F\equiv(mB_{\theta}/r+kB_{z})/B$. Assuming an exponential time
dependence results in a coupled eigenvalue problem with a second-order radial
derivative operator acting on $\psi$. As shown in figure 4, numerical
solutions of these equations confirm that they are sufficient to reproduce the
features of the hose instability observed in the solutions of equations 1-9.

Figure 4: Growth rates and frequencies of the hose mode $m=1$ computed with
the full MHD, the model without bulk plasma flow velocity (eqs. 12-13), and
with $\mathbf{b}=\nabla\psi\times\hat{\mathbf{z}}$ (eqs. 14-15).

The agreement between these three calculations confirms that the flow of the
background plasma is unimportant and that the polarization of the magnetic
field is predominantly transverse. Moreover, the form of equation 14 is
equivalent to the field equation derived by Rosenbluth and Weinberg to
describe the resistive hose modeRosenbluth (1960)Weinberg (1967).

The hose mode in these calculations does not extend beyond a given radius.
Figure 5 shows radial profiles of the hose mode. It is clear that there is
some radius $\tilde{r}$ outside which the profiles of the solution for
$\lambda$ and $\psi$ drop to zero. The value of $\tilde{r}$ appears to be
determined by the solution to $\omega+c_{r}F(\tilde{r})=0$.

Figure 5: Radial profiles of the most unstable eigenfunction of equations 14
and 15 for $(m=1,k=-0.1)$ and $\eta=10^{-2}$.

Although there is no rational surface in the tearing sense of $F(r)=0$, the
equation for $\lambda$ becomes singular in the rotating frame of the eigenmode
at the radial location where $\omega+c_{r}F(r)=0$.

We have shown that resistive hose-type modes associated with a RE beam
interacting with a thermal background plasma grow faster than the typical
resistive MHD tearing mode at large resistivity. The scaling of the hose mode
is linear with the resistivity in the low resistivity regime, and the
frequency of the mode is much larger than the growth rate. Since the hose mode
is dominant when the background plasma resistivity is high, it may appear in
post-TQ tokamak plasmas.

For the RE beam profile considered here, the hose mode grows faster than the
tearing mode number at a normalized resistivity of $\eta/aV_{A}\approx 2\times
10^{-4}$. Ignoring for now whether or not this current density profile is
representative of the experiment, we may use the plasma parameters from a
DIII-D discharge reported in Reference Lvovskiy  _et al._ (2023) to estimate
the temperature where this transition occurs. Using $a=0.5,V_{A}\approx
6.6\times 10^{6}$, and $Z_{eff}=3$, the Spitzer formula produces a temperature
estimate of $T_{e}\approx 1.5$eV. This is in the regime of temperatures
reported for the post-TQ plasmas observed in DIII-D. Even in cases where the
hose mode is subdominant to resistive MHD instabilities, it can introduce an
$m=1$ component into the fluctuation spectrum that would be otherwise
unexpected from resistive MHD alone in cases where the minimum $q$-value is
larger than unity.

Future work remains to find cases in which an analytic dispersion relation may
be found to compare with numerical results. We also intend to address the
nonlinear dynamics of the hose mode and more thoroughly characterize its
influence on RE beam formation and termination in tokamaks.

###### Acknowledgements.

We wish to thank Dr. Chang Liu for providing data on the tearing mode results
with runaways. This work is supported by the US DOE through grant DE-
SC00180001.

## References

  * Breizman  _et al._ (2019) B. N. Breizman, P. Aleynikov, E. M. Hollmann, and M. Lehnen, Nuclear Fusion 59 (2019), 10.1088/1741-4326/ab1822.
  * Lehnen  _et al._ (2015) M. Lehnen, K. Aleynikova, P. Aleynikov, D. Campbell, P. Drewelow, N. Eidietis, Y. Gasparyan, R. Granetz, Y. Gribov, N. Hartmann, E. Hollmann, V. Izzo, S. Jachmich, S.-H. Kim, M. Kofçan, H. Koslowski, D. Kovalenko, U. Kruezi, A. Loarte, S. Maruyama, G. Matthews, P. Parks, G. Pautasso, R. Pitts, C. Reux, V. Riccardo, R. Roccella, J. Snipes, A. Thornton, and P. de Vries, Journal of Nuclear Materials 463, 39 (2015).
  * Avinash and Kaw (1988) Avinash and P. Kaw, Nuclear Fusion 28 (1988), 10.1088/0029-5515/28/1/009.
  * Liu  _et al._ (2020) C. Liu, C. Zhao, S. C. Jardin, A. Bhattacharjee, D. P. Brennan, and N. M. Ferraro, Physics of Plasmas 27 (2020), 10.1063/5.0018559.
  * Zhao  _et al._ (2020) C. Zhao, C. Liu, S. C. Jardin, and N. M. Ferraro, Nuclear Fusion 60 (2020), 10.1088/1741-4326/ab96f4.
  * Helander  _et al._ (2007) P. Helander, D. Grasso, R. J. Hastie, and A. Perona, Physics of Plasmas 14 (2007), 10.1063/1.2817016.
  * Rosenbluth (1960) M. N. Rosenbluth, Physics of Fluids 3, 932 (1960).
  * Weinberg (1964) S. Weinberg, Journal of Mathematical Physics 5, 1371 (1964).
  * Weinberg (1967) S. Weinberg, Journal of Mathematical Physics 8, 614 (1967).
  * Fernster  _et al._ (1995) R. F. Fernster, S. P. Slinker, M. Lampe, and R. F. Hubbard, Physics of Plasmas 2, 4338 (1995).
  * Uhm and Lampe (1980) H. S. Uhm and M. Lampe, Physics of Fluids 23, 1574 (1980).
  * Lvovskiy  _et al._ (2023) A. Lvovskiy, C. Paz-Soldan, N. W. Eidietis, A. D. Molin, G. H. DeGrandchamp, E. M. Hollmann, J. B. Lestz, C. Liu, M. Nocente, D. Shiraki, and X. D. Du, Nuclear Fusion 63 (2023), 10.1088/1741-4326/acb725.
  * Sovinec  _et al._ (2004) C. Sovinec, A. Glasser, T. Gianakon, D. Barnes, R. Nebel, S. Kruger, D. Schnack, S. Plimpton, A. Tarditi, and M. Chu, Journal of Computational Physics , 355 (2004).
  * Sovinec (2016) C. R. Sovinec, Journal of Computational Physics 319, 61 (2016).
  * Hollmann  _et al._ (2013) E. M. Hollmann, M. E. Austin, J. A. Boedo, N. H. Brooks, N. Commaux, N. W. Eidietis, D. A. Humphreys, V. A. Izzo, A. N. James, T. C. Jernigan, A. Loarte, J. Martin-Solis, R. A. Moyer, J. M. Muñoz-Burgos, P. B. Parks, D. L. Rudakov, E. J. Strait, C. Tsui, M. A. V. Zeeland, J. C. Wesley, and J. H. Yu, Nuclear Fusion 53 (2013), 10.1088/0029-5515/53/8/083004.


# Don’t do it: Safer Reinforcement Learning With Rule-based Guidance

Ekaterina Nikonova,
Cheng Xue, Jochen Renz

###### Abstract

During training, reinforcement learning systems interact with the world
without considering the safety of their actions. When deployed into the real
world, such systems can be dangerous and cause harm to their surroundings.
Often, dangerous situations can be mitigated by defining a set of rules that
the system should not violate under any conditions. For example, in robot
navigation, one safety rule would be to avoid colliding with surrounding
objects and people. In this work, we define safety rules in terms of the
relationships between the agent and objects and use them to prevent
reinforcement learning systems from performing potentially harmful actions. We
propose a new safe epsilon-greedy algorithm that uses safety rules to override
agents’ actions if they are considered to be unsafe. In our experiments, we
show that a safe epsilon-greedy policy significantly increases the safety of
the agent during training, improves the learning efficiency resulting in much
faster convergence, and achieves better performance than the base model.

## Introduction

In reinforcement learning, the agent tries to maximize the long-term return
and often explores the world in a random manner without considering the safety
of its actions. While safety is usually less important in simulated
environments, in the real world, safe interaction is often a rather expected
and necessary quality of the agent. For example, for a delivery robot it is
important to not collide with any surrounding obstacles damaging itself or the
package; for a self-driving car making safe decisions is crucial for the
safety of the driver and others. In AI safety, this particular issue is known
as safe exploration [6] and has been widely explored in reinforcement learning
[1, 2, 3, 4, 21, 22, 5].

Multiple approaches have been taken by the researchers to provide safer
reinforcement learning. For example, transforming optimization criteria [7, 8,
9] or providing external guidance [12, 10, 11, 13].

Often, safe exploration and interaction can be ensured by stating safety rules
for the particular domain. For example, consider the earlier mentioned real-
world task for the delivery robot to navigate the streets without crashing
into surrounding objects and people. In this task, the most important safety
rule to follow is to avoid performing the action that would result in an
immediate crash. This rule can be easily defined using spatial relationships
between the robot and other objects, for example using qualitative spatial
representation [19, 18]. In words, we can define a simple rule for this task
as ”if the object is in close proximity to the agent and the intended action
would result in a collision - do not perform such action”.

Figure 1: (Left) General framework to validate the safety of the agent’s
decision using safety component. (Right) An example of framework
implementation in Freeway. Here predicting action ”up” would result in
violating rule 1 and such a decision would be overwritten by safer action
”down”.

The idea of stating rules to prevent agents from doing something unsafe goes a
long way back to Issac Asimov’s ideas in his ”Three Laws of Robotics” [24]. In
some cases, such basic safety rules are either known and can be injected as
prior knowledge to the agent, or can be learned by the agent as it explores
the world.

In reinforcement learning, researchers have explored the idea of using prior
knowledge to address safe exploration problem [1, 21, 22]. One possible way to
use prior knowledge is by using imitation learning [23]. For example, the
agent could imitate human demonstrations to avoid performing unsafe actions
[25]. Another approach is to include the safety constraints directly in the
policy optimization [26]. Saunders proposed to entirely override the agent’s
actions if it was seen to be unsafe [27]. We find Saunders’s idea to be the
most similar to ours, however, they use human oversight to override agent
decisions and learn intervention policy from human demonstrations. While shown
to be effective and safe, their approach requires a human to be in the loop
which is both time-consuming and often unrealistic. In our work, we extend
this idea and eliminate the need of constant human oversight. Instead, we use
either human-defined or machine-learned safety rules to override agents’
actions if they are unsafe. To the best of our knowledge, we have not found a
method that uses symbolic safety rules directly to prevent an agent from
performing unsafe actions.

In this work, we focus on the domains that require agents to act safely to do
a task well. We propose a general framework to integrate the safety component
together with reinforcement learning and provide one possible implementation
of that framework using symbolic rules (Figure 1). In our experiments, we use
human knowledge about the internal dynamics of the environment to manually
construct symbolic safety rules. To show that our framework can work with
different safety components, we provide a short experiment to learn the rules
as the agent explores the world. We propose a new method to incorporate safety
rules into the reinforcement learning algorithm. Specifically, we define a
safe-$\epsilon$-greedy policy that uses symbolic rules to validate the safety
of actions and overrides them with safer actions if necessary. In all our
experiments we use the same underlying deep learning model with identical
hyperparameters to ensure a fair evaluation of our method. We demonstrate that
restricting the agent in its actions improves the overall performance,
increases the safety of the agent, and significantly improves the efficiency
of the agent achieving super-human or state-of-the-art performance faster. In
our experiments, we test our method in five domains: a toy-domain Crossroad,
and four OpenAI Gym Atari 2600 environments: Freeway, Frostbite, Ms-Pacman,
and Space Invaders [14].

## Background

In this work we consider a deterministic Markov Decision Process (MDP)
$M=(S,A,T,r,\gamma)$, where $S$ is the state space, $A$ is the action space,
$T:S\times A\rightarrow S$ the transition function, $r:S\times A\rightarrow r$
is a reward function, and $y\in[0,1)$ the discount factor.

A policy $\pi:S\rightarrow A$ determines which action to take in each state.
Typically, the goal of reinforcement learning is to find a policy that
maximizes the expected discounted reward and is therefore considered to be
optimal.

The Q-function
$Q^{\pi}(s,a)=\mathbb{E}^{\pi}[\sum^{\infty}_{t=0}\gamma^{t}r_{t}|s_{0}=s,a_{0}=a]$
measures the performance of the agent assuming it starts in a state $s$, takes
action $a$ and follows the policy $\pi$ afterwards.

The Value-function $V^{\pi}(s)=\mathbb{E}^{a\sim\pi(s)}[Q^{\pi}(s,a)]$
measures the overall value of the state. Same as with policy, those functions
can be optimal: $Q^{*}(s,a)=max_{\pi}Q^{\pi}(s,a)$ and
$V^{*}(s)=max_{\pi}V^{\pi}(s)$. Finally, the optimal policy can be retrieved
from $Q^{*}$ as follows: $\pi^{*}(s)=argmax_{a}Q^{*}(s,a)$.

In deep reinforcement learning, $Q$-function can be approximated using a
nonlinear function approximator such as a neural network $Q(s,a,\theta_{i})$,
where $\theta_{i}$ are the weights of the Q-network at the i-th iteration.
However, when using a nonlinear function approximator together with
reinforcement learning, it can become unstable or even diverge due to the
following problems: a) the correlation in the sequence of observations, b)
correlations between Q values and target values $r_{t}+\gamma
max_{a}Q(s_{t},a)$ and c) having a policy that is extremely sensitive to
changes of Q value.

A deep Q-network (DQN) addresses the first problem by using experience replay.
Experience replay is implemented by storing and later randomly sampling the
observations experienced by the agent. This technique removes the correlation
between the sequences of the observations by randomizing the collected data.
We define the experience as $e_{t}=(s_{t},a_{t},r_{t+1},s_{t+1})$, and
experience set as $M=\\{e_{1},…,e_{t}\\}$.

In order to address the second problem, the notion of target network was
introduced which is then used to calculate the loss:
$L_{i}(\theta_{i})=\mathbb{E}_{(s,a,r,s^{{}^{\prime}})\sim U(M)}[(r+\gamma
max_{a^{{}^{\prime}}}Q(s^{{}^{\prime}},a^{{}^{\prime}},\theta_{i}^{-})-Q(s,a,\theta_{i}))^{2}]$,
where $i$ is the iteration, $\gamma$ is a discount factor, $\theta_{i}$ are
weights of so-called online Q-network and $\theta_{i}^{-}$ are weights of a
target network or so-called offline Q-network. The target network is called
offline since its weights are only updated every $C$ steps with a copy of
online network weights, while the online network is updating every iteration
$i$.

## Guidance Using Symbolic Rules

### Definitions

#### Qualitative Spatial Representation

In order to construct a symbolic representation of the states, we extract
qualitative spatial relationships (QSR) between the objects and the agent. In
particular we use cone-shaped directional representation [29] and a
qualitative distance representation [30]. Figure 2 shows an example of a
combination of directional and distance representations. We use this
combination to construct safety rules and relational representation of the
state.

Figure 2: Qualitative directional and distance representation used for the
safety rules.

#### Safety Rules

One obvious approach would be to define safety in terms of the state
transitions and next states [2]. However, such a definition would require an
agent to have a model which is not always possible. Instead, we look at the
relationships between the objects to determine if the action is safe or not in
the given state. We define a safety rule to be a conjunction of n-ary
relationships between objects and action which if satisfied would compromise
safety: $safetyRule=(r_{1}(o_{1},...o_{n})\wedge...\wedge
r_{m}(o_{1},...,o_{n})\wedge\lnot action(a))$, where $r_{i}$ is a QSR
relationship, $o_{j}$ is an object and $action(a)$ is the action. For example,
a safety rule would prevent agent from performing an action that could lead to
the collision with the near-by objects (i.e. $close(agent,o)\wedge
N(agent,o)\wedge\lnot action(up)$). We define safety rules to be a collection
of such rules. Figure 1 (Right) demonstrates symbolic safety rules for
Freeway.

#### Is the action safe?

We call an action $a_{t}$ to be safe in state $s_{t}$, if performing that
action would not violate any safety rule. To validate that action is safe, we
extract symbolic relationships from the state $s_{t}$ as
$s^{symb}_{t}=\\{r_{1}(o_{1},...o_{n}),...,r_{m}(o_{1},...,o_{n})\\}$, where
$r_{i}$ is a QSR relationship and $o_{j}$ is an object in the state $s_{t}$.
For example, we can extract spatial relationships between the nearby objects
and the agent (i.e. $close(agent,car_{1})\wedge N(agent,car_{1})$). We then
conjugate it with the symbolic representation $a^{symb}_{t}$ of action $a_{t}$
(i.e. $action(a_{t})$). Finally, we define $isActionSafe(s_{t},a_{t})$ =

$\begin{cases}1,&\text{if }(s^{symb}_{t}\wedge\lnot a^{symb}_{t})\cap
safetyRules=\emptyset\\\ 0,&\text{otherwise}\end{cases}$ (1)

Where $s^{symb}_{t}$ are the relationships extracted from $s_{t}$ and
$a^{symb}_{t}$ is symbolic representation of $a_{t}$. We note that in this
work, we only extract spatial relationships between objects and agents that
fall within a certain region around the agent (Figure 1 (Right)). Consider
Figure 1 that demonstrates a framework implementation in a Freeway. Here the
safety rules prevent the action that would result in a collision with a car.
In this example the car is immediately above the agent (i.e.
$s^{symb}=(close(agent,car_{1})\wedge N(agent,car_{1}))$). If the DQN predicts
action ”up” (i.e $close(agent,car_{1})\wedge N(agent,car_{1})\wedge
action(up)$), rule 1 would be violated, and such action would be overridden
with a safe action.

Algorithm 1 Guided Exploration

Input: $Q$, $s_{t}$
Output: $a$

1: n $\sim$ $\mathcal{U}_{[0,1]}$

2: if n <$\epsilon$ then

3: a $\leftarrow$ selectRandomSafeAction($s_{t}$)

4: else

5: a $\leftarrow$ $argmax_{a}Q(s_{t},a)$

6: end if

7: return a

#### Random safe action

To select a random safe action in the state $s_{t}$, we define
$selectRandomSafeAction(s_{t})=$

$\begin{cases}a\in A_{safe},&\text{if }A_{safe}\neq\emptyset\\\ a\in
A,&\text{otherwise}\end{cases}$ (2)

In Equation 2, action is sampled from $A_{safe}$ or $A$ uniformly and
$A_{safe}=\\{a^{i}_{t}|isActionSafe(s_{t},a^{i}_{t})=1,a^{i}_{t}\in A\\}$ is a
collection of all safe actions. The actions that would result in violating any
safety rule are called to be unsafe and should be avoided by the agent
whenever it is possible.

### Safe $\epsilon$-greedy policy

In this section, we propose a method to inject the safety rules as part of the
modified, safer version of the $\epsilon$-greedy algorithm. Safe
$\epsilon$-greedy algorithm directly prevents the agent from performing unsafe
action by overriding its decision with a safe action (if possible). We
hypothesize that restricting the agent in its actions, will not only make
learning safer but likely will reduce the size of the search space and
therefore improve the performance.

In each time-step, safe $\epsilon$-greedy leverages safety rules to determine
the safety of agents’ actions. We provide two variants of the safe
$\epsilon$-greedy algorithm: Guided exploration and Full-guidance. First
validates and overrides agents’ actions during the exploration phase only,
while second during the whole training process. We describe the algorithms
further in the section.

Recall the $\epsilon$-greedy policy [16]:

$\pi(s_{t})=\begin{cases}a\in$A$,&\text{if
}n\textless\epsilon,n\in\mathcal{U}_{[0,1]}\\\
argmax_{a}Q(s_{t},a),&\text{otherwise}\end{cases}$ (3)

By definition, the epsilon-greedy policy (Equation 3) selects a random action
during the exploration phase, which, without any safety check, could result in
the agent damaging itself or others. Instead, we propose to select a random
safe action.

Figure 3: Testing domains used in this work. From left to right: Crossroad,
Freeway, Frostbite, Ms-Pacman, and Space Invaders. Algorithm 2 Full Guidance

Input: $Q$, $s_{t}$
Output: $a$

1: n $\sim$ $\mathcal{U}_{[0,1]}$

2: if n <$\epsilon$ then

3: a $\leftarrow$ selectRandomSafeAction($s_{t}$)

4: else

5: a $\leftarrow$ $argmax_{a}Q(s_{t},a)$

6: if not isActionSafe($s_{t}$, $a$) then

7: a $\leftarrow$ selectRandomSafeAction($s_{t}$)

8: end if

9: end if

10: return a

Algorithm 1 demonstrates how safety rules can then be embedded into the
$\epsilon$-greedy policy resulting in its safer version. Here instead of
selecting a random action, the algorithm selects random safe action. We note
that safe action will be added to the experience replay $M$ (or Q-table as in
tabular RL) to aid further training (feedback in Figure 1).

While Algorithm 1 will ensure the agent is behaving safely during the
exploration, it does not guarantee that once the exploration phase is over the
agent will continue to behave safely. For example, if the agent encounters a
previously unseen state, it could perform an unsafe action. Therefore, to
ensure that the agent behaves safely during the entire training process, we
must continue to supervise its actions.

Algorithm 2 presents an approach that provides the agent with guidance during
the whole training process. Same as Algorithm 1, it will select safe actions
during the exploration phase. However, in addition to it, when an action is
selected using a learned $Q$-function, the algorithm will check whether this
action is safe, and if not, it will randomly select a safe action. We
hypothesize that such strict supervision of the actions performed by the agent
will increase the safety of the agent. Ultimately, as the result of this
algorithm, many unsafe actions will not be explored thus reducing the size of
the search space. One could note that theoretically doing something unsafe
could be a necessary step to achieve maximum performance, however here we
prioritize safety over reward.

## Experiments

#### Testing Domains

In this work, we use five testing domains: Crossroad, and four OpenAI Gym
Atari 2600 domains: Freeway, Frostbite, Ms-Pacman, and Space Invaders (Figure
3).

Crossroad is an initial testing domain inspired by the Atari 2600 Freeway
game. It has 7 cars moving horizontally either left or right and discrete
action space with the ability to move left, right, down, and up or stay in
place. The goal of the game is to cross the road without colliding with a car.
Colliding with a car will result in the immediate end of the game. The agent
is given the reward of 1 if the game is won, -1 if it has collided with a car,
and 0 otherwise.

Atari 2600 Freeway, is a game similar to Crossroad with cars moving left or
right horizontally. However, in comparison to Crossroad, agents can only move
up and down or stay in place. The goal of the game is to cross the road the
maximum possible number of times before the timer runs out. Here colliding
with a car will not result in the end of the game. The agent is given the
reward of 1 if it crosses all roads and 0 in any other case. In Atari 2600
Frostbite, the agent needs to jump onto the platforms to build an igloo while
avoiding falling into the water. In Atari 2600 Ms-Pacman agent needs to
collect all dots while avoiding the ghosts (however, ghosts can be eaten in
some cases for extra reward). Lastly, in Atari 2600 Space Invaders the agent
needs to shoot enemies while dodging the bullets.

We chose these five domains as to win and master those games, the agent needs
to act safely. For example, in Crossroad and Freeway agent needs to avoid the
cars, in Ms-Pacman it should avoid ghosts, in Frostbite it should not jump
into the water and in Space Invaders, it should avoid the bullets. All of this
can be reformulated to real-world tasks such as avoiding obstacles and not
falling off from high places. In addition to those five domains, our method
can be easily applied to any other domain with well-defined safety rules.

#### Agents and Settings

We evaluate safe $\epsilon$-greedy policy in tabular and deep reinforcement
learning settings. As representation learning is not a focus of this work, we
use symbolic state representation for the Crossroad and RAM for Atari 2600. In
order to extract symbolic representation for the Atari games, we use the AARI
package [15]. We note, however, that this part can be easily replaced with any
state-of-the-art image processor to extract the location of the objects. We
test two variants of the safe $\epsilon$-greedy policy using Guided-
exploration and Full-guidance algorithms.

Tabular Q-learning allows us to measure the impact of the method without
complications and instability of deep reinforcement learning. For tabular
Q-learning, we test three agents: vanilla Q-learning, Q-learning + Guided
Exploration, and Q-learning + Full-guidance.

Figure 4: Average rewards of DQN agents in Crossroad with 0-1 reward.

For deep reinforcement learning experiments, we use Deep Q-network as a base
model due to its relative simplicity and widespread in the research community.
However, our method can be used with any other off-policy algorithm as a
substitute for $\epsilon$-greedy policy. We test three agents: vanilla DQN,
DQN + Guided Exploration (DQN+GE), and DQN + Full-guidance (DQN+FG). We note
that we use the same hyperparameters and model architecture for all three of
our tested agents. The only difference between the three agents is in action-
selection algorithm they use during training (either Equation 3, or Algorithms
1 and 2).

We compare the performance of the agents in Atari games to the performance of
humans, random agent [28], and image-based DQN (image-DQN) that was trained
for 50 million frames [17]. We include the results of image-based DQN to show
the difference between the RAM and image state representation.

#### Safety Rules for Testing Domains

As mentioned earlier, safety rules can be defined in terms of the
relationships between the objects and actions. In this work, we define
symbolic safety rules for each testing domain manually. However, we note that
manually defining the rules is both cumbersome and time-consuming [20]. To
address this, we hypothesize that those rules can be learned autonomously.
Later in the paper (in ”Learning Rules Experiment” section), we present a
short experiment showing that: a) rules can be learned and b) using learned
rules will work with our method as well.

In both Crossroad and Freeway the agent should not collide with the cars. We
define safety rules in terms of relationships between the cars and the agent
(Figure 1 (Right)). For example, to determine whether the action up is safe we
check whether there is a car directly above the agent. If not then it is safe
to move up. Similar rules apply to all other directions and actions.

In Frostbite the agent should avoid falling into the water. Similarly, safety
rules can be defined in terms of the relationships between platforms and the
agent. For example, to determine whether the action ”down” is safe, we check
whether there is a platform beneath the agent.

Figure 5: Average rewards of Q-learning agents and learning rule agent in
Crossroad.

In Ms-Pacman, the rule is quite similar to Crossroad and Freeway - to not
collide with the ghosts, with the exception that sometimes ghosts can be eaten
to receive bonus points. As previously this rule can be defined in terms of
relationships between the game objects.

Lastly, in Space Invaders, the agent needs to avoid the bullets, which can be
defined in terms of the relationships between the bullets and the player. When
the bullet is close to the player, the player needs to move to dodge it. We
provide full safety rules for all domains in the Appendix.

#### Difference to the Reward Shaping

One logical question that one may have is whether constraining the agent in
its action will lead to similar performance as providing the negative rewards
for such actions. One obvious downside of this approach is that a negative
reward will not completely prevent the agent from performing unsafe actions.
The agent would still need to perform an action to learn that it is unsafe.
Often, in the real world, this would be too much of a risk. We hypothesize,
that since the agent would still need to explore the larger search space,
likely it will not outperform the safe $\epsilon$-greedy agents.

Figure 4 provides the learning curves of DQN agents in Crossroad. Originally,
as for example in Figure 8, Crossroad provides the reward of -1 if the agent
collides with the car, 1 if it crosses all roads, and 0 in all other cases.
However, we modified the reward function for this experiment to give a reward
of 1 for crossing all roads and 0 otherwise. DQN+Neg.Reward is the only agent
that was receiving a negative reward when colliding with a car, thus emulating
the scenario where human shapes the reward function for each domain. In this
experiment, the agent with a negative reward did not outperform the safe
$\epsilon$-greedy agents nor vanilla DQN with 0-1 reward.

#### Learning Rules Experiment

One obvious limitation of using human-defined safety rules is that it takes
time to define and implement them [20]. However, we hypothesise that such
rules can be learned by agents automatically and conduct a short experiment to
test this hypothesis. We demonstrate that learned rules can boost efficiency
similarly to manually defined rules. We note, that learning rules is beyond
the focus of this work and leave the full-scale experiment to future work.

Figure 6: The scaled number of times each agent has either died or collided
relatively to the base model such as DQN or Q-learning. The lesser the number
the safer training is (meaning fewer collisions or deaths).

For this experiment we use supervised learning as safety component in our
framework (Figure 1 (Left)). We train logistic regression on the collection of
symbolic states and classify them to ”safe” and ”unsafe” categories. In each
state $s_{t}$ we use the model to get the next state $s_{t+1}$ following the
action $a_{t}$ and use classifier to determine whether the $s_{t+1}$ is safe.
We provide full details on the classifier and the approach in the Appendix.
Figure 5 shows the learning curves of the four agents: vanilla Q-learning,
Q-learning with learned rule and Guided exploration, and two agents that use
human-defined rules: Q-learning + Guided exploration, and Q-learning + Full-
exploration. We conclude that using a learned rule is virtually almost
identical to human-defined rules. This concludes that rules can be learned and
used in safe $\epsilon$-greedy policy, we leave the full-scale experiment with
other domains for future work.

## Results and Discussions

In this subsection, we present the results of the experiments in five domains.
We discuss how a safe $\epsilon$-greedy policy increases the safety of the
agent during training, improves the learning efficiency, and achieves better
performance than the base models.

#### Increased Safety

Safe $\epsilon$-greedy leverages the safety rules to override agents’ actions
if they are considered to be unsafe. We hypothesize that this would decrease
the number of collisions or deaths during training and potentially help the
agent to learn a safer policy. Figure 6 shows the relative number of times
each agent has died (or collided) during training. In all domains safe
$\epsilon$-greedy agents demonstrated a noticeable reduction in the number of
deaths. Based on those results, we hypothesize that in the real world our
method would reduce the potential cost of training as the agent would damage
itself fewer times. We note that in some domains the decrease in the deaths
was not as significant as in other domains. We hypothesise that this is due to
the short horizon of the safety rules. In some domains (i.e. Ms-Pacman), agent
needs to plan its actions several steps ahead to avoid states where death is
inevitable. One possible way to deal with it would be to use a more
sophisticated safety component to validate the safety of the agents’ actions.

Figure 7: The number of episodes it took each agent to achieve the maximum
performance of the vanilla DQN or Q-learning. The lesser the number the more
efficient the agent was (meaning faster convergence).

#### Improved Efficiency

As we mentioned previously, we hypothesize that by preventing the agent from
exploring unsafe actions, we could decrease the size of the search space thus
increasing convergence speed. Figure 7 demonstrates how many episodes it took
the safe $\epsilon$-greedy agents to achieve maximum performance of the base
model. It demonstrates that for all domains the efficiency has improved
significantly, achieving good performance faster than the base model. Figure 8
demonstrates average rewards achieved by agents in all domains. Safe
$\epsilon$-greedy policy provided the agents with the initial jump-start in
the performance and, in some cases, increased the convergence speed. Here we
note that for Crossroad, vanilla DQN failed to learn after 25,000 episodes,
whereas both DQN+GE and DQN+FG have converged to the optimal performance. The
initial jump-start resulted in both agents surpassing vanilla DQN performance
in the first episode.

Figure 8: Average rewards of the agents in five testing domains.  Agent | Crossroad | Freeway | Frostbite | Ms-Pacman | Space Invaders
---|---|---|---|---|---
DQN | -0.8 | 27.3 | 358.0 | 2780.0 | 1155.0
DQN+GE | 1.0 | 28.9 | 490.0 | 3558.0 | 1470.0
DQN+FG | 1.0 | 30.5 | 470.0 | 3820.0 | 1670.0
Random | -1.0 | 0 | 65.2 | 307.3 | 148.0
image-DQN | - | 30.3 | 328.3 | 2311.0 | 1449.7
Human | 1.0 | 29.6 | 4334.7 | 6951.6 | 1668.7
Table 1: The maximum reward achieved after training for 25,000 episodes (<20M
frames). Here image-DQN are the results of image-based DQN reported in [17].

#### Better Performance

In addition to bringing safety and efficiency, a safe $\epsilon$-greedy policy
also improves the performance of the agents, resulting in achieving a higher
maximum reward. Table 1 demonstrates the maximum rewards achieved by the
agents during training in all testing domains. We compare the achieved maximum
performance of the agents to humans and a random agent. We also provide the
results of the image-based DQN [17] to compare the performance impact of two
different representations. In all domains safe $\epsilon$-greedy agents
achieved better performance than the base models. We hypothesize that this is
likely due to the two main reasons: increased efficiency and initial jump-
start in the performance, and chosen testing domains. The first reason is
rather self-explanatory, as better efficiency would result in achieving better
performance faster, and an initial jump-start would provide a better starting
point for the learning. Earlier, we mentioned that we selected the domains
that require an agent to behave safely in order to achieve good performance.
We hypothesize, that in comparison to some other domains, here safety
correlates with the performance. Since the agent explores only safe states, it
finds better policy faster.

## Conclusions and Limitations

In this work, we proposed to use domain-specific symbolic safety rules to
prevent agents from performing unsafe actions. We proposed a new safe
$\epsilon$-greedy algorithm that injects symbolic rules into the training
process of the reinforcement learning agent. We presented two versions of this
algorithm that either completely restricts the agent in its actions during the
entire training process, or only during the exploration of the world. We
demonstrated how the safety rules can be defined for five different domains.
In our experiments, we showed that preventing an agent from performing unsafe
actions can significantly increase the efficiency of the agent, boost overall
performance and make the training process much safer. We provided an
additional experiment to show how eliminating unsafe actions completely
results in much better learning efficiency rather than simply giving a
negative reward. Finally, we addressed a known limitation and showed that
safety rules can be learned and integrated into the training using our method.

One obvious limitation of using symbolic safety rules is that in more complex
domains rules can be much more complicated and hard to define. However, as
mentioned in the ”Learning Rules Experiment” section, we hypothesize that in
such complicated domains the rules can be learned by the agent using known
rule-learning techniques. Since this goes beyond this paper, we leave it for
future work.

Another limitation is that our method does not completely eliminate unsafe
actions (Figure 6). This is mainly due to the short horizon of defined
symbolic rules. The rules that we use are only able to model one step ahead.
This leaves a number of situations where inevitably agent has to perform some
unsafe action. For example, in Freeway, there could be situations where the
agent does not have enough time or any way to avoid colliding with a car. To
avoid such scenarios more sophisticated and complex rules are needed. We
hypothesize that to avoid unsafe actions completely a model may be needed,
which essentially leads to a realm of planning. Another way could be to use a
neurosymbolic approach. However, this discussion goes beyond this paper and we
leave this to future work.

In conclusion, safe-$\epsilon$ greedy policy provides a flexible way to inject
safety rules into the learning process to improve the safety and efficiency of
the reinforcement learning agents.

## References

  * [1] García, J. & Fernández, F. A comprehensive survey on safe reinforcement learning. J. Mach. Learn. Res.. 16 pp. 1437-1480 (2015)
  * [2] Hans, A., Schneegaß, D., Schäfer, A. & Udluft, S. Safe exploration for reinforcement learning. ESANN. (2008)
  * [3] Thomas, G., Luo, Y. & Ma, T. Safe Reinforcement Learning by Imagining the Near Future. NeurIPS. (2021)
  * [4] Kobelrausch, M. & Jantsch, A. Collision-Free Deep Reinforcement Learning for Mobile Robots using Crash-Prevention Policy. 2021 7th International Conference On Control, Automation And Robotics (ICCAR). pp. 52-59 (2021)
  * [5] Gu, S., Yang, L., Du, Y., Chen, G., Walter, F., Wang, J., Yang, Y. & Knoll, A. A Review of Safe Reinforcement Learning: Methods, Theory and Applications. ArXiv. abs/2205.10330 (2022)
  * [6] Amodei, D., Olah, C., Steinhardt, J., Christiano, P., Schulman, J. & Mané, D. Concrete Problems in AI Safety. ArXiv. abs/1606.06565 (2016)
  * [7] Geibel, P. & Wysotzki, F. Risk-Sensitive Reinforcement Learning Applied to Control under Constraints. J. Artif. Intell. Res.. 24 pp. 81-108 (2005)
  * [8] Howard, R. & Matheson, J. Risk-Sensitive Markov Decision Processes. Management Science. 18, 356-369 (1972),
  * [9] Lütjens, B., Everett, M. & How, J. Safe Reinforcement Learning With Model Uncertainty Estimates. 2019 International Conference On Robotics And Automation (ICRA). pp. 8662-8668 (2019)
  * [10] Chow, Y., Nachum, O., Duéñez-Guzmán, E. & Ghavamzadeh, M. A Lyapunov-based Approach to Safe Reinforcement Learning. NeurIPS. (2018)
  * [11] Clouse, J. & Utgoff, P. A Teaching Method for Reinforcement Learning. ML. (1992)
  * [12] García, J. & Fernández, F. Safe Exploration of State and Action Spaces in Reinforcement Learning. ArXiv. abs/1402.0560 (2012)
  * [13] Geramifard, A., Redding, J. & How, J. Intelligent Cooperative Control Architecture: A Framework for Performance Improvement Using Safe Learning. Journal Of Intelligent and Robotic Systems. 72 pp. 83-103 (2013)
  * [14] Brockman, G., Cheung, V., Pettersson, L., Schneider, J., Schulman, J., Tang, J. & Zaremba, W. OpenAI Gym. ArXiv. abs/1606.01540 (2016)
  * [15] Anand, A., Racah, E., Ozair, S., Bengio, Y., Côté, M. & Hjelm, R. Unsupervised State Representation Learning in Atari. ArXiv. abs/1906.08226 (2019)
  * [16] Sutton, R. & Barto, A. Reinforcement Learning: An Introduction. (A Bradford Book,2018)
  * [17] Mnih, V., Kavukcuoglu, K., Silver, D., Rusu, A., Veness, J., Bellemare, M., Graves, A., Riedmiller, M., Fidjeland, A., Ostrovski, G., Petersen, S., Beattie, C., Sadik, A., Antonoglou, I., King, H., Kumaran, D., Wierstra, D., Legg, S. & Hassabis, D. Human-level control through deep reinforcement learning. Nature. 518, 529-533 (2015,2), http://dx.doi.org/10.1038/nature14236
  * [18] Clementini, E., Felice, P. & Hernández, D. Qualitative Representation of Positional Information. Artif. Intell.. 95 pp. 317-356 (1997)
  * [19] Cohn, A. & Renz, J. Qualitative Spatial Representation and Reasoning. Handbook Of Knowledge Representation. (2008)
  * [20] Weld, D. & Etzioni, O. The First Law of Robotics (A Call to Arms). AAAI. (1994)
  * [21] Pecka, M. & Svoboda, T. Safe Exploration Techniques for Reinforcement Learning - An Overview. MESAS. (2014)
  * [22] Leike, J., Martic, M., Krakovna, V., Ortega, P., Everitt, T., Lefrancq, A., Orseau, L. & Legg, S. AI Safety Gridworlds. ArXiv. abs/1711.09883 (2017)
  * [23] Santara, A., Naik, A., Ravindran, B., Das, D., Mudigere, D., Avancha, S. & Kaul, B. RAIL: Risk-Averse Imitation Learning. AAMAS. (2018)
  * [24] Asimov, I. I, Robot. (Conn: Fawcett Publications,1950,6)
  * [25] Ghavamzadeh, M., Petrik, M. & Chow, Y. Safe Policy Improvement by Minimizing Robust Baseline Regret. NIPS. (2016)
  * [26] Achiam, J., Held, D., Tamar, A. & Abbeel, P. Constrained Policy Optimization. ICML. (2017)
  * [27] Saunders, W., Sastry, G., Stuhlmüller, A. & Evans, O. Trial without Error: Towards Safe Reinforcement Learning via Human Intervention. ArXiv. abs/1707.05173 (2018)
  * [28] Lazaridis, A., Fachantidis, A. & Vlahavas, I. Deep Reinforcement Learning: A State-of-the-Art Walkthrough. J. Artif. Intell. Res.. 69 pp. 1421-1471 (2020)
  * [29] Renz, J. & Mitra, D. Qualitative Direction Calculi with Arbitrary Granularity. PRICAI. (2004)
  * [30] Frank, A. Qualitative spatial reasoning about distances and directions in geographic space. J. Vis. Lang. Comput.. 3 pp. 343-371 (1992)


# Solid State Neuroscience: Spiking Neural Networks as Time Matter

Marcelo Rozenberg  Université Paris-Saclay, CNRS Laboratoire de Physique des
Solides, 91405, Orsay, France

###### Abstract

We aim at building a bridge between to a priori disconnected fields:
Neuroscience and Material Science. We construct an analogy based on
identifying spikes events in time with the positions of particles of matter.
We show that one may think of the dynamical states of spiking neurons and
spiking neural networks as time-matter. Namely, a structure of spike-events in
time having analogue properties to that of ordinary matter. We can define for
neural systems notions equivalent to the equations of state, phase diagrams
and their phase transitions. For instance, the familiar Ideal Gas Law relation
(P$v$ = constant) emerges as analogue of the Ideal Integrate and Fire neuron
model relation ($I_{in}$ISI = constant). We define the neural analogue of the
spatial structure correlation function, that can characterize spiking states
with temporal long-range order, such as regular tonic spiking. We also define
the “neuro-compressibility” response function in analogy to the lattice
compressibility. We show that similarly to the case of ordinary matter, the
anomalous behavior of the neuro-compressibility is a precursor effect that
signals the onset of changes in spiking states. We propose that the notion of
neuro-compressibility may open the way to develop novel medical tools for the
early diagnose of diseases. It may allow to predict impending anomalous neural
states, such as Parkinson’s tremors, epileptic seizures, electric
cardiopathies, and perhaps may even serve as a predictor of the likelihood of
regaining consciousness.

##  I Introduction

The understanding of the mind is, arguably, the most mysterious scientific
frontier. The ability of the mind to understand itself is puzzling.
Nevertheless, it seems increasingly possible and within our reach.
Neuroscience and Artificial Intelligence are making great progress in that
regard, however, following very different paths and driven by very different
motivations. In the first, the focus is to address fundamental questions of
biology, while in the second, it is to develop brain inspired computational
systems for practical applications for modern life. Evidently, there is also
large overlap between them both.

The basic units that conform the physical support of the mind, namely the
brain and the associated neuronal systems, are neurons. These are cells with
electrical activity that interact via electric spikes, called action
potentials [1]. In animals, neurons conform networks of a wide range of
complexity. Ranging from a few hundred units in jelly fish, to a hundred
billion in humans. A fundamental question to answer is how and why nature has
adopted this electric signaling system. Its main functions are multifold: to
sense and monitor the environment, then to produce behavior and decision
making, and finally to drive the required motor actions to assure the survival
of living beings. Neuroscience has already provided a good understanding of
the electric behavior of individual neurons [2, 3]. A major milestone was the
explanation by Hodgkin and Huxley of the physiological mechanism for the
generation of the action potential [4]. At the other end, that of neural
networks with large number of neurons, important contributions come from
Artificial Intelligence. For instance, significant progress was made in the
80s following the pioneering work of Hopfield [5]. More recently, this area
received a renewed boost of activity, enabled by the combination of new
learning algorithms for Deep Convolutional Neural Networks [6] with the
numerical power of modern computers [7]. However, the networks adopted in
Artificial Intelligence overwhelmingly describe the neurons’ activity by their
firing rate and not by individual spikes. These are conceptually different, a
spike is a discrete event, while the spiking rate is a continuous variable.
Hence, modeling neurons in terms of the latter does not directly address the
question posed above, namely, the why and how of Nature using discrete spikes.

Here we propose to look at this problem under a different light, which to my
knowledge has not been discussed before. Hopefully, this may bring new
insights and perhaps help to develop our intuition for the challenging problem
of understanding the mechanism of spiking neural networks. We shall postulate
an analogy between matter states, as in organized spatial patterns of
particles (or atoms or molecules) and that of neuronal states, as organized
patterns of spikes in time. Since we are attempting to build a bridge between
disconnected disciplines, we shall keep our presentation pedagogical.

Ultimately, as with any definition, our analogy will be of value if it turns
to be a useful besides the intrinsic academic interest. With this in mind, we
shall later discuss an exciting perspective that the present approach may
perhaps open. Namely, to develop novel screening tools that could allow to
detect an enhanced risk to develop neural diseases that involve anomalous
spiking states, such as Parkinson’s, epilepsy, cardiopathies, unconsciousness.

##  II The analogy

As mentioned above, here we postulate an analogy between spatial matter
states, such as solid, liquid and gas, and the dynamical states of spiking
neurons and neural networks. More specifically, an analogy between the
organization of matter in space and that of spiking events in time, which we
call temporal-matter states.

To motivate this, in Fig.1 we show the familiar phase diagram of water as
function of pressure and temperature ($p$,$T$). Next to it we show another
phase diagram, that of an electronic bursting spiking neuron, which we
introduce recently [8]. In the diagram, we observe various phases that
correspond to qualitatively different states: tonic spiking (TS), fast spiking
(FS), and two types of bursting (IB1, IB2). The phase diagram is obtained as a
function of two parameters, the excitatory current $I$ and a circuit time
constant $\tau_{s}$ (the circuit is shown in Fig.4).

Figure 1: Left: Phase diagram of states of ordinary matter (water) showing the
solid, liquid and gas phases. $P$ and $T$ determine the specific volume $v$ of
the state, which is the inverse of the density $n$. Right: Measured phase
diagram of temporal-matter states of an electronic neuron model. The system is
a two-compartment bursting neuron, showing tonic spiking (TS), fast spiking
(FS) and two different intrinsic bursting states (IB1 and IB2). The input
current $I$ and the circuit time-constant $\tau_{s}$ determine the inter-spike
interval of the state, which corresponds to the inverse of the frequency.

To characterize the spiking states in the phase diagram, we need to consider
the nature of the discrete spiking events in the time domain. For instance,
tonic spiking is characterized by a sequence of spikes that occur at equally
spaced time-intervals. In Neuroscience, the time between two consecutive spike
events is called the inter-spike interval (ISI). At the transition from the TS
to FS phase, one observes a sudden decrease of the ISI (i.e. a jump in the
spiking frequency) as function of the parameter $\tau_{s}$.

The ISI(t) characterizes the organization of the spikes in the time domain and
it indicates the “time-distance” between spikes. If a sequence of spike events
at times $t_{i}$ are indicated by the function

$s(t)=\sum_{i}\delta(t-t_{i})$ (1)

Then, it may be tempting to establish the analogy by simply identifying the
time-position of spike events with the positions of particles in a matter
state. The matter state is characterized by its particle density function,

$n(x)=\sum_{i}\delta(x-x_{i})$ (2)

where $x_{i}$ denotes the position of the $i^{th}$ particle. For simplicity,
we assume point particles and one dimensional space. Similarly, in Eq.1 we
have assume ideal spikes represented by $\delta$-functions, while in reality
the action potentials have a typical duration of $\sim 1ms$ From the analogy,
the simple tonic spiking case with equispaced spikes (i.e. constant ISI) would
correspond to a perfect crystal of equally spaced particles (atoms). Thus, the
familiar ISI of neuroscience would be the analogue of the familiar lattice
constant $a$ (or specific volume $v$) for material science.

If one applies a pressure $P$ to matter, in general one observes the decrease
of $v$, such as in the case of the Ideal Gas law $P$$v$ = constant. On the
other hand, in neuroscience, it is well known that the ISI can be reduced by
increasing the excitatory input current. Hence, one may be tempted to extend
our analogy to associate the $P$ with the input $I$ in neural systems. We can
make this more precise by introducing the simplest theoretical model of an
ideal spiking neuron, the Integrate and Fire (IF) model [2].

In a very schematic view, as shown in Fig.2, a biological neuron is composed
of three main parts: the dendrites, the soma and the axon. The neuron is
excited through the input of electric signals arriving to the dendrites, which
are called the synaptic currents. This input is integrated in the cell’s body,
which leads to the increase in its electric potential with respect to its
resting state. By means of intense excitation, the neuron eventually reaches a
threshold potential value. At that point a dramatic event takes place, an
electric spike is initiated, it propagates down the axon and eventually is
communicated to the dendrites of a downstream neuron. This phenomenon is
called the emission of an action potential, which was first described by
Hodgkin and Huxley [4]. This qualitative description can be represented by the
simple (leaky) IF model [2, 3]

$\frac{du}{dt}=-\frac{1}{\tau}u(t)+I;\ \ if\ u\geq u_{th}\ then\ spike\ and\
u=u_{rest}$ (3)

where $u(t)$ represents the potential of the soma, $u_{th}$ is the threshold
potential, $u_{rest}$ is the resting potential value and $I$ is the input
(synaptic) current. The time constant $\tau$ is a characteristic relaxation
time of the neuron that represents the leakage of charge out of the soma. If
the leakage is negligible, $\tau\to\infty$, and the integration is perfect, so
one has an ideal IF model. The electric circuit representation of the model is
straight forward. The soma is represented by a capacitor $C$ that accumulates
the charge of the input current, and the leakage by a resistor $R$ in
parallel. The threshold voltage can be represented by a switch that closes
yielding the emission of the spike, which is the fast discharge of the charge
accumulated in $C$, as a delta function of current.

Figure 2: (a) Schematic biological neuron. (b) Electric circuit of the IF
model (for the case $u_{res}$=0). The membrane of the cell’s body (soma) is
represented by the capacitor $C_{m}$, which accumulates charge of the input
current. If the resistor $R\to\infty$ then the current integration is perfect,
however for a finite value there is a “leaky” integration and the model is
called LIF, for leaky-integrate-and-fire.

The limit of zero leakage, i.e. the ideal Integrate and Fire is trivially
solved. The potential $u$ due to the integrated charge in $C$ during an
interval $t$ is $u(t)=Q/C=(I/C)t$. Thus, the spike fire time $t_{f}$ is given
by the condition $u(t_{f})=u_{th}$. Which leads to $It_{f}=u_{th}C$, or $I{\rm
ISI}=constant$.

Hence, we may extend our analogy by noting that the equation of state of the
Ideal Gas has the same form as that of an ideal IF neuron, namely,

$Pv{\rm=constant}\longleftrightarrow I{\rm ISI=constant}$ (4)

where

$P\longleftrightarrow I;\ \ \ v\longleftrightarrow{\rm ISI}$ (5)

As the equation of state of a real gas departs from the ideal case, biological
neurons and neuron models’ equation of state, will depart from the ideal IF
“neuronal equation of state” above. Interestingly, the notion of neuronal
equation of state should not appear so strange. Indeed, it is nothing else
than the familiar concept in neuroscience of the neuron’s activation function,
namely, the firing rate as a function of the excitatory input current,
$f=f(I)$. For instance, some popular models are: the rectified linear units,
or ReLU, where $f(I)=max[0,I]$; the sigmoid activation $f(I)=1/(1+e^{-I})$;
etc. These are examples of mathematical neuron models, however, we may also
include here physical neuron models, namely, models that are defined by an
electronic circuit. In physical neuron models the equation of state $f(I)$ can
be measured, as in real gases. We should mention that while the relation $Pv$=
constant is familiar from the Ideal Gas law, it is in general valid for any
liquid or solid in the linear regime.

Recently, we introduced a minimal model of a physical spiking neuron, which is
achieved by exploiting the memristive properties of an “old” conventional
electronic component, the thyristor. The model is minimal because it provides
a physical realization of the basic leaky-integrate-and-fire (LIF) neuron
model by associating exactly one component to each one of the three functions:
a capacitor to integrate, a resistor to leak and the thyristor to fire.
Qualitatively, the thyristor acts as the switch in the circuit of Fig.2.

Figure 3: (a) Electric circuit that defines the physical neuron model based on
the memristive properties of the thyristor. The function of the thyristor is
that of a voltage controlled switch, as shown in Fig.2 above. The small “load”
resistor $R_{L}$ is to transform the output spike current into a voltage
action potential signal for measurement convenience. (b) The “neuron equation
of state” or activation function $f(I)$. Near the threshold, the equation of
state follows closely the functional form of a LIF model (red line),
$-1/[log(1-I_{c}/I)]$ where $I_{c}$ is the minimal activation current. The
blue line corresponds to the Ideal IF behavior $I\propto$ 1/ISI$\ =f$.

In Fig.3 we show the circuit that defines physical neuron model just
described, where we identify the role of each of the three components to the
functions of the LIF. We call this artificial neuron model the Memristive
Spiking Neuron (MSN), which is implicitly defined by its electronic circuit.
In the right panel of the figure we show the experimental neuron equation of
state, which is nothing other than the activation function, as noted above. We
can observe that near the excitation threshold the equation of state is well
represented by the functional form of the activation function of the LIF
mathematical model (red fitting line in Fig.3). At intermediate input
currents, the behavior approaches the Ideal IF neuron, whose equation of state
is $f\propto I$ (see 4), since $f=1/ISI$

The same methodology allows us to consider more complex neuron models, which
are also defined by their respective circuit implementation. For instance, we
may consider the case of bursting neurons. From theoretical neuroscience we
know that a requirement to obtain bursting behavior is by adding a second
dynamical variable, besides the potential of the cell body $u(t)$, in Eq.4 and
represented by the capacitor in the circuits of Fig.2 and 3. Thus, to do this
we add a second $RC$ pair to our basic MSN circuit. A simple option is to add
a capacitor $C_{L}$ in parallel to the (small) output resistor $R_{L}$,
introducing a new time constant $\tau_{L}=R_{L}C_{L}$. The resulting circuit
is shown in Fig.4a. In panel (b) we show the dynamical behavior that it now
produces. We observe four qualitatively different spiking types: simple tonic
spiking, fast spiking, and two bursting modes. These four spiking types are
realized in the respective regions of the phase diagram of Fig.1 presented
before.

Figure 4: (a) Electric circuit that defines the physical model of a bursting
neuron based on the MSN model and adding a second time constant
$\tau_{L}=R_{L}C_{L}$. (b) The various dynamical behaviors produced by the
circuit. From top to bottom: tonic spiking (TS), fast spiking (FS) and two
bursting types (IB1, IB2), which correspond to the four phases of the phase
diagram shown in Fig.1.

As done before for the basic MSN, we may also obtain the equation of state of
the Memristor Bursting Neuron (MBN) model. In the present case, we can
consider that the time constant $\tau_{L}$ plays the role of a third
parameter, similarly as the temperature in the case of matter systems. Hence,
in Fig.5 we show the curve $f_{\tau_{L}}(I)$ measured at a fixed $\tau_{L}$,
indicated by the vertical purple line that crosses three phases. We observe
the jumps in the frequency as the current drives the system from one phase to
the other. This is reminiscent of changes in density when the pressure drives
phase transitions at a fixed $T$ in the phase diagram of water in Fig.1.

It is interesting to mention that a biologically realistic theoretical model
of bursting neurons introduced by Pinsky and Rinzel (PR) shows a qualitatively
similar behavior [9]. The PR is an example of a two-compartment model, where
both the soma and the dendrites are described. We may notice that in the MBN
model the $R_{L}C_{L}$ block (see Fig.4) can be considered as a second
compartment, which is connected to the output of the first. In the right panel
of Fig.5, we reproduce the activation function (i.e. the neuronal equation of
state) of the PR model. It is interesting to observe that in the simpler
limits of only one compartment (soma alone and dendrite alone) the behavior of
the PR is qualitatively the same as that of our basic MSN, which is also a
single compartment model (see bottom curve of Fig.5a). More importantly, for
the relevant case of two compartments (i.e. finite $g_{c}$) we observe a
sudden changes in the firing rate as a function of excitatory input current,
also in qualitative agreement with the MBN. In fact, both PR and MBN traverse
the same sequence as the excitatory current is increased: initially quiescent
below a critical current, then bursting, and finally jump in firing rate to
the fast spiking mode. Hence, the phase transitions are abrupt, through a
steep or a discontinuous change in the activation function. As we shall
discuss in the next section this feature may have interesting consequences.
Moreover, and we shall show that the $f(I)$ anomalies may be considered as the
counterparts of certain phenomena occurring in phase transitions in matter
systems.

Figure 5: (a) Activation function $f(I)$ in semilog scale. The top curve
corresponds to the MBN model along the right purple line at $\tau_{L}$=0.3ms
indicated in the phase diagram in the inset. Following the definition of
Pinsky and Rinzel, the frequency is defined as the inverse of the period
between spike trains (bursts). The bottom $f(I)$ curve is that of the MSN
model discussed before, which corresponds to the limit $\tau_{L}\to 0$ of the
bursting neuron model (indicated with the left purple line in the inset). (b)
Activation function of the Pinsky-Rinzel model reproduced from [9].

##  III Correlation and Response Functions

Correlation and response functions are useful concepts in material science,
where they serve to characterize different states of matter. For instance, a
regularity in the arrangement of positions of atoms is revealed by Bragg peaks
in the x-ray spectra, which are maxima of the structure factor. In real space,
the regularity is revealed by the pair correlation function

$g({x})=\frac{\int n(x+x^{\prime})n(x^{\prime})dx^{\prime}}{\int
n^{2}(x^{\prime})dx^{\prime}}$ (6)

where $n(x)$ indicates the particle density (such as electrons, atoms,
molecules, etc) at position $x$, and where we consider one dimension for
simplicity. In the case of a crystalline order, $g(x)$ shows structure with
peaks. For a simple arrangement of particles along one dimension with a
lattice constant $a$, the peaks will be at $a$, $2a$, $3a$, … In contrast, for
a disordered state, such as gas or liquid, the $g(x)$ is mostly featureless.
The study of $g(x)$ is routinely done in condensed matter physics for the
study of phase transitions (see, for example, [10]).

In our analogy, spiking systems are thought of as temporal-matter states, so
it is natural to explore the behavior of the correlation-function analogue to
$g(x)$. Since the position of particles correspond to position of spikes, by
analogy we can define the neural correlation function $g_{n}(t)$ as

$g_{n}(t)=\frac{\int s(t+t^{\prime})s(t^{\prime})dt^{\prime}}{\int
s^{2}(t^{\prime})dt^{\prime}}$ (7)

where $s(t)$ indicates a given spike trace. The function $g_{n}(t)$ can
characterize different spiking states of a neuron. In Fig.6 we provide a
concrete example, which is realized in the basic MSN model described before
(see Fig.3).

Figure 6: (a) Measured spike traces $s(t)$ of the MSN model at two values of
the input current: a “solid” state measured at $I$ = 45.3$\mu$A (top) and a
“melted” state at $I$ = 36.4$\mu$A. The states are indicated with the green
and blue dots in the activation function of the neuron reproduced in the
inset. (b) The neural correlation function $g_{n}(t)$ computed for the
respective traces shown in the left side panels (only a small portion of the
measured traces is shown).

We can observe the two qualitatively different tonic spiking behaviors in the
two left panels of the figure. They correspond to two constant current inputs
(36.4$\mu$A and 45.3$\mu$A). In the first case, at higher input current (top
panel), the spiking is perfectly regular. In contrast for a smaller current
close to the threshold, the trace changes dramatically. The interval between
spikes become very irregular. By analogy between spike and particles, we can
think of the first case as that of a solid and the second as of the melting of
the solid state. This qualitative description can be made more precise by the
correlation function $g_{n}(t)$, shown on the right side panels of Fig.6.

The top panel shows a succession of delta functions at equally spaced at
times, multiples of the constant inter-spike interval, $t_{k}$=$k$ISI. This
indicates the long range order in time. It tells that given the presence of a
spike at time $t=0$ we have a high probability ( 1) of finding another spike
event at times $t_{k}$ ($k$=1, 2, 3, …). In contrast, the $g_{n}(t)$ shown in
the bottom panel is featureless, showing a small approximate constant value.
This indicates total lack of order as the presence of a spike a time $t=0$
does not permit to predict the presence of ulterior spikes events. The
emission of spikes is random as the positions of particles in a liquid or a
gas.

The peaks of the $g_{n}(t)$ are very narrow, delta-function-like, because the
spikes are very narrow with respect to the duration of the ISI. In a solid,
the atoms have a size that is smaller but of the same order as the lattice
spacing, so instead of narrow deltas one observes broad peaks in the $g(x)$
[10].

One may understand quite intuitively the physical origin of this dramatic
changes in the time-structure. The key point is to realize that the “melted”
state occurs in a regime where the activation function $f(I)$ is very steep,
at the onset of neuron excitability, i.e. near the threshold. Therefore, small
variations of the input current will reflect on significant variations in the
ISI. This observation motivates the following important insight.

This enhanced sensitivity to current fluctuations is due to the large slope
$df/dI$ in the activation function. Then, what is this feature related to, if
one follows the analogy back to the matter systems? We recall that ISI plays
the role of the lattice spacing, or specific volume (see Eqs.4 and 5), then
$f$=1/ISI corresponds the particle density $n(x)=1/v(x)$. On the other hand,
since the input current $I$ is like the pressure, then it follows that the
slope $df/dI$ corresponds to $dn/dP$. This last quantity is closely related to
the compressibility of matter systems

$\beta=\frac{1}{n}\ (\frac{dn}{dP})$ (8)

which is the inverse of the bulk modulus. We can therefore follow the analogy
and introduce the concept of “neuro-compressibility”,

$\beta_{n}=\frac{1}{f}\ (\frac{df}{dI})$ (9)

It is important to mention that this quantity may be measured using
experimental methods such as dynamic clamp, where a controlled synaptic
current can be injected into a neuron while its activity is monitored [11].
Moreover, this definition may turn out to have important consequences, as we
discuss next.

Anomalies in the compressibility of materials are precursor signatures of
structural phase transitions. A sudden increase in the compressibility of a
solid indicates the “softening” of a vibrational mode (a phonon mode), which
leads to a change in the structure, or possibly a phase change. Then, the
question is, what would the analogue phenomenon be for a neuronal system? For
a single neuron, the enhancement of $\beta_{n}$ would indicate the proximity
to a qualitative change in the spiking mode of a neuron, i.e. a “bifurcation”
in its dynamics. This can in fact be seen in the panels of Fig.5. There, we
observe that all the changes in the spiking modes for both, the MSN and for
the MBN models, occur at current values where there are enhancements in
$df/dI$ or jumps. Most notably, this is not only a feature of our artificial
neuron circuits, but also can be clearly seen in the biologically realistic
Pinsky-Rinzel model activation functions that we reproduced in Fig.5 [9]. The
enhancements seen in the PR data occur at the onset of change from quiescent
to spiking and also in the change from burst to tonic spiking (black circles
and black triangles), in very good qualitative agreement to our electronic
neuron model,

We may then speculate on an important implication of our observations. It
would be interesting to explore if neuro-compressibility anomalies are also
found across the boundaries of qualitatively different states in neuron
networks. If that is the case, an intriguing and exciting possibility would be
to investigate if anomalies in $\beta_{n}$ are also detected (by small current
stimulation) in animal models of epilepsy and Parkinson’s disease. If this
were the case, then one may envision a pathway to a novel diagnose tool for
early detection or a risk predictor of mental diseases associated to abnormal
spike patterns in humans. In even further speculation, one may also search for
anomalies at the onset of regaining, or loosing, consciousness, which is
another challenging frontier of research [12].

##  IV Bursting spikes as an ananlogue of formation of clusters defects

Here we describe another interesting connection between common a phenomena in
spiking neurons and in material science. We shall show that missing spikes in
a trace of a fast spiking state, can be thought of as the analogue of missing
atoms, i.e, defects, in a crystal structure.

From Fig.7 we observe that the proliferation of missing spikes in a fast
spiking state is a route to generate bursting behavior. This is illustrated in
the sequence of traces shown in the figure, which were obtained for a step-
wise decreasing input current to the MBN. The thick purple arrow in the
vertical path followed in phase diagram (from blue to green to grey). In the
top trace we indicate with small purple arrows the missing arrows, showing
that one may understand the onset of bursting as the result of skipping spike
events, which are initially few (i.e. dilute). As the current intensity is
further reduced, the missing spikes become more numerous (i.e dense) and
occurring in clusters of inactivity, which give rise to the stuttering mode
bursts [13]. In our analogy, we think of spikes as of atoms in a lattice,
therefore, the initial continuous fast spiking state is like a perfect
crystal. The missing spikes then play the analogue role as vacancy defects,
i.e. lacking atoms. Moreover, the missing spikes are the result of decreasing
current, which in the analogy represents pressure. It is then interesting to
observe that in thin-film deposition, which is a topic in material science,
the partial pressure of oxygen $P$(O2) is a relevant parameter for the quality
of the growth of crystalline oxides. Moreover, it is well known that reducing
the $P$(O2) induces the creation of oxygen vacancy defects in the crystal
structure [14, 15], which often cluster together forming dislocations [16],
This is in full qualitative analogy to the spiking traces in the stuttering
bursting mode shown in Fig.7.

We would like to emphasize that the path of phase transformations, i.e. the
evolution from fast spiking, to bursting, to quiescent does not seem to be
just a peculiarity of our MBN circuit model. In the lower panel of the Fig.7
we illustrate the striking resemblance of the traces of the MBN with those
measured in bursting neuron of rats [17]. Quite remarkably, the experimental
traces were obtained by solely changing the intensity of the excitatory DC
current.

Figure 7: Top left: Measured spike traces of the MBN as a function of
decreasing current in discrete steps (purple). The thin purple arrows indicate
the missing spikes. Bottom left: Experimental trace of pre-Bötzinger bursting
neurons from rats. The data is obtained by changing excitatory input current
in discrete steps. Adapted from [17]. The circuit parameters may be easily
adjusted to fit the experimental data [8]. Right: Phase diagram of the MBN
where the purple arrow indicates the evolution from the fast spiking (blue) to
the bursting type 1 phase (green) by decreasing current at constant
$\tau_{L}$.

##  V Neural Networks

We now consider one important final aspect in our analogy that may eventually
bring new light to the issue of how to think about inter-neuron coupling. So
far we have considered essentially individual neurons, but we may ask what
would it be to extend the analogy to multi-neuron systems, i.e. to neural
networks.

As a first glimpse into this question, we shall consider the simplest network
case, namely, just two neurons that are mutually excitatory or inhibitory. We
focus first in a system of two identical neurons, each excited with equal
input currents. The currents are above threshold, so the neurons are active
and theirs spikes are transformed via conductances into mutually injected
synaptic currents that are positive for the excitatory case and negative for
the inhibitory one. This is schematically shown in Fig.8.

We shall see that our analogy takes an interesting twist, as the dynamical
states of the two-neuron system can be considered as an analogue of a complex
crystal, i.e. a crystal with two atoms in the unit cell. Moreover, the
coupling between neurons is mediated by synaptic currents, which can be
excitatory or inhibitory. Since in our analogy current plays the role of
pressure, then the synaptic currents should also play such a role. More
precisely, an excitatory synaptic current should correspond to a repulsive
inter-particle interactions (positive pressure), and an inhibitory synaptic
current should be an attractive interaction (negative pressure). We shall
consider the two cases, which we study using realistic electronic circuit
simulations (LTspice).

Figure 8: Top: two 1D space lattices (red and blue) with repulsive interaction
and two neurons with inhibitory interactions. The neurons are prepared to
start spiking synchronously and spontaneously adopt an alternating spiking
pattern. Bottom: two 1D lattices with attractive interactions and two neuron
network with excitatory interaction. They are prepared to alternate but
spontaneously begin to spike in unison. Spike traces were produced by
realistic simulations with LTSpice.

We first consider two spiking neurons with mutual inhibition. When one neuron
spikes it inhibits the other, and vise-versa, so they try to avoid spiking at
unison since we are synaptic current are instantaneous. After a transient
period, they expectedly find a stable dynamical state where they alternate to
emit spikes as shown in the top panel of Fig.8. Using our analogy, this
spiking pattern corresponds two a perfect molecular crystal where the unit
cell has an A-B atom pair (or a basis).

In the bottom panel of the figure we consider the second case, that of
mutually excitatory neurons. Again after a transient time, the system adopts a
periodic spiking pattern. However, in contrast to the previous case, now both
neurons fire at unison. This is also intuitive, the excitatory synaptic
current of emitted by a neuron that spiked promotes the spiking of the other
one, and vise versa. So, naturally, the both spike at the same time, which is
nothing else that excitatory synapses promote synchrony in neural networks
[18]. Following our analogy, spiking at unison corresponds to a “dimerization”
of the lattice. Namely, that the distance between the A-B pair atoms is
reduced as due to an attractive interaction between the A and B atomic
species.

These two cases are consistent with our analogy, where current is interpreted
as pressure. Indeed, the volume of the unit cell in the inhibitory case is
large, as expected for positive effective pressure between A and B, while in
the excitatory case the volume is fully collapsed to zero, as expected for a
negative pressure within the unit cell.

It is an interesting perspective for future work to consider increasingly
complex networks of several neurons (motifs). The periodic states that emerge
constitute spiking sequences, which are of great relevance for automatic motor
behavior. By virtue of the analogy that we introduced in the present work,
those periodic sequences should correspond to a variety of molecular crystals.
It would be exciting to explore if new intuition for Neuroscience could be
brought from from those traditional areas of Condensed Matter Physics and
Chemistry [19].

##  VI Conclusion

In this work we introduced the idea that the dynamical states of neural
networks may be though of as realizations of “time matter” states.

We started from the notion that a trace of spiking events of a neuron can be
analogue to a snapshot of particles or atoms arranged in space. We then went
on to explore and show that the analogy may be pushed far beyond that literary
statement, and may provide new intuition in the challenging problem of
understanding and designing spiking neural networks.

We identified analogue roles of basic quantities of Physics with those of
Neuroscience, such as pressure and volume with input currents and inter-spike
intervals. We then logically build on this assumption to show connections
between correlation and response functions in both fields. Perhaps most
significant was the finding that a neuro-compressibility can be defined, with
possible far reaching consequences, including medical ones, that may be
experimentally tested.

An exciting new road of discovery may open ahead.

##  VII Acknowledgments

We acknowledge support from the French ANR “MoMA” project ANR-19-CE30-0020.

## References

  * Bear  _et al._ [2016] M. Bear, B. Connors, and M. Paradiso,  _Neuroscience: Exploring the Brain_ (Wolters Kluwer, Philadelphia, 4Th ed. international ed., 2016).
  * Gerstner  _et al._ [2014] W. Gerstner, W. Kistler, R. Naud, and L. Paninski,  _Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition_ (Cambridge University Press, 2014).
  * Ermentrout and Terman [2010] G. Ermentrout and D. Terman,  _Mathematical Foundations of Neuroscience_, Interdisciplinary Applied Mathematics (Springer New York, 2010).
  * Hodgkin and Huxley [1952] A. L. Hodgkin and A. F. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve, The Journal of Physiology 117, 500 (1952).
  * Hopfield [1982] J. J. Hopfield, Neural networks and physical systems with emergent collective computational abilities, Proceedings of the national academy of sciences 79, 2554 (1982).
  * LeCun  _et al._ [2015] Y. LeCun, Y. Bengio, and G. Hinton, Deep learning, Nature 521, 436 (2015).
  * Merolla  _et al._ [2014] P. A. Merolla, J. V. Arthur, R. Alvarez-Icaza, A. S. Cassidy, J. Sawada, F. Akopyan, B. L. Jackson, N. Imam, C. Guo, Y. Nakamura, _et al._ , A million spiking-neuron integrated circuit with a scalable communication network and interface, Science 345, 668 (2014).
  * Wu  _et al._ [2023] J. Wu, O. Schneegans, K. Wang, P. Stoliar, and M. Rozenberg, submitted (2023).
  * Pinsky and Rinzel [1994] P. F. Pinsky and J. Rinzel, Intrinsic and network rhythmogenesis in a reduced traub model for ca3 neurons, Journal of Computational Neuroscience 1, 39 (1994).
  * Debela  _et al._ [2014] T. T. Debela, X. D. Wang, Q. P. Cao, Y. H. Lu, D. X. Zhang, H.-J. Fecht, H. Tanaka, and J. Z. Jiang, Nucleation driven by orientational order in supercooled niobium as seen via ab initio molecular dynamics, Phys. Rev. B 89, 104205 (2014).
  * Goaillard and Marder [2006] J.-M. Goaillard and E. Marder, Dynamic clamp analyses of cardiac, endocrine, and neural function, Physiology 21, 197 (2006), pMID: 16714478\.
  * Storm  _et al._ [2017] J. F. Storm, M. Boly, A. G. Casali, M. Massimini, U. Olcese, C. M. A. Pennartz, and M. Wilke, Consciousness regained: Disentangling mechanisms, brain systems, and behavioral responses, The Journal of Neuroscience 37, 10882 (2017).
  * Golomb  _et al._ [2007] D. Golomb, K. Donner, L. Shacham, D. Shlosberg, Y. Amitai, and D. Hansel, Mechanisms of firing patterns in fast-spiking cortical interneurons, PLOS Computational Biology 3, e156 (2007).
  * Ghenzi  _et al._ [2015] N. Ghenzi, M. J. Rozenberg, R. Llopis, P. Levy, L. E. Hueso, and P. Stoliar, Tuning the resistive switching properties of tio2-x films, Applied Physics Letters 106, 123509 (2015).
  * Liu  _et al._ [2009] X. Liu, K. Zhou, L. Wang, B. Wang, and Y. Li, Oxygen vacancy clusters promoting reducibility and activity of ceria nanorods, Journal of the American Chemical Society 131, 3140 (2009).
  * Szot  _et al._ [2006] K. Szot, W. Speier, G. Bihlmayer, and R. Waser, Switching the electrical resistance of individual dislocations in single-crystalline srtio3, Nature Materials 5, 312 (2006).
  * Butera  _et al._ [1999] R. J. J. Butera, J. Rinzel, and J. C. Smith, Models of respiratory rhythm generation in the pre-bötzinger complex. i. bursting pacemaker neurons., J Neurophysiol 82, 382 (1999).
  * Hansel  _et al._ [1995] D. Hansel, G. Mato, and C. Meunier, Synchrony in Excitatory Neural Networks, Neural Computation 7, 307 (1995).
  * Corpinot and Bučar [2019] M. K. Corpinot and D.-K. Bučar, A practical guide to the design of molecular crystals, Crystal Growth & Design 19, 1426 (2019).