A speed-accurate self-sustaining head direction cell path integration model without recurrent excitation

ABSTRACT

The head direction (HD) system signals HD in an allocentric frame of reference. The system is able to update firing based on internally derived information about self-motion, a process known as path integration. Of particular interest is how path integration might maintain concordance between true HD and internally represented HD. Here we present a self-sustaining two-layer model, capable of self-organizing, which produces extremely accurate path integration. The implications of this work for future investigations of HD system path integration are discussed.

KEYWORDS:

• Spatial cognition
• head direction cells
• path integration
• continuous attractor neural networks
• attractor dynamics

Author Summary

Head direction (HD) cells are neurons that represent the current HD of an animal with respect to the environment. Whilst not geomagnetic, they represent an internal compass. HD cells can anchor to the environment on the basis of visual information, but can also use vestibular and other self-motion signals to update their firing as the animal’s head rotates. This process is known as path integration.

Using computer simulation, we investigate a new model of the rat HD system which learns to perform highly accurate path integration (when the internal representation of HD matches the current direction of the head). Our model also attempts to explain two further related questions: why the operation of the HD system seems to be dependent on vestibular input, and how HD cells, without direct excitatory projections to one another, can maintain their firing in the absence of visual input.

We propose a system in which cells responding to conjunctions of HD and angular head velocity can be used to both sustain HD cell firing and update this firing at the correct speed during head turns. We also demonstrate how the required synaptic connectivity between neurons can develop during an initial period of visually-guided learning.

Introduction

Head direction (HD) cells respond to the animal’s HD in the yaw plane (Taube et al. 1990a). Their responses reflect a combination of externally derived (allothetic) and internally derived (idiothetic) information (Knight et al. 2014; Page et al. 2014). Although they are anchored to visual cues (Taube et al. 1990b), HD cells are able to sustain and update their firing based solely on idiothetic information, a process known as path integration (Mittelstaedt and Mittelstaedt 1980; Etienne and Jeffery 2004).

One important question relating to the path integration process is how concordance is maintained between the current HD representation and true HD. That is, how does the HD system learn to perform path integration at close to the correct speed? Although the system may accumulate error, which is likely to be reset by visual landmarks, path integration is quite accurate (Goodridge et al. 1998). As it currently stands, this is an underinvestigated question. One theory of accurate path integration (Walters and Stringer 2010; Walters et al. 2013) involves associations being formed between successively active HD cells via connectivity containing some natural time delay (e.g., axonal conduction delay). In effect, this allows for prediction of where an HD activity packet will be after a fixed time delay, given a constant rotational velocity and a known starting location. HD cells project to a postsynaptic HD cell with an offset in HD space determined by time delay and rotational velocity. However, computer simulations exploring this hypothesized mechanism have demonstrated at best $81\mathrm{\%}$ accuracy (Walters 2011; Walters et al. 2013).

The investigation into the sources of simulated path integration inaccuracy (Page et al. 2015) reveals two major factors. First is the presence of within-layer recurrent collateral excitatory connectivity between HD cells. Such connectivity can only learn a single rotational speed of any arbitrary value (including 0). Recurrent connectivity will, therefore, interfere with the system’s attempts to update HD cell firing at anything other than this rotational speed. Second is the time between when a neuron receives stimulation and when it begins firing, known as rise time ($t_R$). In a system in which path integration speed is based on a specific conduction delay, the rise time represents an increase in this delay. Specifically, if a postsynaptic HD cell is supposed to fire some time interval $\mathrm{\Delta }t$ after presynaptic firing, then it will, in fact, begin firing at $\mathrm{\Delta }t+t_R$. In this manner, rise time introduces error in path integration.

This paper addresses the first of these issues. Following the attractor hypothesis of Skaggs and colleagues (Skaggs et al. 1995), most computational models have been based on continuous attractor neural networks (CANNs). However, such models typically incorporate within-layer excitatory recurrent collaterals of the sort which will negatively affect path integration accuracy. In addition, there is currently no evidence for such connectivity in the HD system in vivo. There have been several HD system models without such recurrent connectivity, although they appear to have been motivated by this latter physiological reason rather than a consideration of path integration accuracy (Rubin et al. 2001; Boucheny et al. 2005; Song and Wang 2005).

The model in this paper extends previous work by providing a self-organizing two-layer HD system model capable of highly accurate path integration. We begin by demonstrating near-perfect path integration accuracy in a pre-wired model with fixed synaptic weights in order to clarify the computational mechanisms underpinning accurate path integration. We then extend this work to a self-organizing system using biologically plausible local associative learning rules to set-up the required synaptic connectivity. The implications of these results on the computational understanding of path integration inaccuracy are discussed.

Materials and methods

Pre-wired model description

The model in this paper draws from the existing literature of HD cell modelling and can be considered an extension of the standard CANN approach. It consists of three layers: one layer of HD cells, one layer of cells responding to angular head velocity (AHV cells), and one layer of cells responding to combinations of HD and AHV (COMB cells). Its architecture is pictured in Figure 1. Within the HD ring, each cell i has a preferred firing direction $x_i$ in the range $0-360{}^{\circ }$. These directions are distributed evenly around the HD ring. Cells are arranged topographically such that adjacent HD cells have adjacent preferred firing directions. This arrangement is for ease of visualization and analysis. It does not reflect HD cell organization in the brain, where no such topography has been observed, nor does it affect model operation. In this manner, a contiguous Gaussian packet of activity in the HD layer represents current head direction. A crucial aspect of this model is that HD cells do not have a direct excitatory influence on one another. Instead, HD cells only project to ($w^1$) and receive from ($w^2$) the COMB cell layer. This HD to COMB and COMB to HD connectivity contains an axonal conduction delay $\mathrm{\Delta }t$ in both directions.

Figure 1. Network architecture of pre-wired model.

The model consists of three layers of neurons: a head direction (HD) layer, a layer of COMB cells comprizing ROT COMB and NOROT COMB cells, and a layer of angular head velocity (AHV) cells which include both ROT and NOROT cells. The layers are connected by $w^1$ connections from HD to COMB cells, $w^2$ connections from COMB to HD cells, $w^3$ connections from ROT cells to ROT-COMB cells, and $w^4$ connections from NOROT cells to NOROT-COMB cells.

Cells in the COMB layer receive connections from the HD and AHV layers and project back onto the HD layer. The AHV layer consists of two kinds of a cell: a ROT cell, which is active during head rotation, and a NOROT cell which is active when the head is not rotating. Accordingly, COMB cells are conceptually split into two subpopulations: those that receive from the ROT cells (ROT-COMB cells) via $w^3$ connectivity and those that receive from the NOROT cells (NOROT-COMB cells) via $w^4$ connectivity. Connectivity from the AHV to COMB layer does not contain any delays. Both NOROT and ROT variants of COMB cells have preferred directions similar to HD cells, by virtue of the connectivity they receive from the HD cell layer, as described below. These two COMB cell variants have separate preferred direction distributions, each evenly distributed from $0-360{}^{\circ }$.

Whilst ROT and NOROT cells are programmed as distinct in this model, it should be kept in mind that this is a simplification of the known physiology: in reality, they originate from the same layer of cells, and do not exist as distinct populations. Indeed, the binary nature of the ROT and NOROT cells simply signals whether or not the head is rotating at the target velocity. This too is a simplification of angular head velocity (AHV) cells in the brain, which have complex firing responses to angular head velocity. The aim with these ROT and NOROT cells is simply to capture the fact that the AHV cell population as a whole will provide information about when the head is rotating at any given AHV, including 0, to the COMB cell layer.

Connectivity between the HD and COMB, as well as that between COMB and HD layers, is based on COMB subpopulations. Individual NOROT-COMB cells receive and project from the same HD cell. These NOROT-COMB cells thus respond to a combination of head direction and a lack of head rotation. ROT-COMB cells receive maximally from an offset region of the HD layer and project back onto a further offset region of the HD layer. These offsets are based on the conduction delay $\mathrm{\Delta }t$ in each direction and the speed of head rotation signalled by the corresponding ROT cells. The operating principles of this connectivity will be fully detailed below.

Individual HD and COMB cells are simulated using a leaky-integrator firing rate model. Consequently, individual spikes are not modelled, and instead, cell firing is measured as an instantaneous average firing rate. ROT and NOROT cells are modelled as binary units, being either on (firing rate of 1) or off (firing rate of 0).

Pre-wired model equations

The activation level h${ }_i$${ }^{\mathrm{H}\mathrm{D}}$ of head direction cell i is given by

(1) ${\tau }^{HD}\frac{d{h}_i^{HD}(t)}{dt}=-h_i^{HD}(t)+e_i(t)-\frac{1}{N^{HD}}\sum _j{\tilde{w}}^{HD}{r}_j^{HD}(t)+\frac{{\varphi }_2}{C^{COMB\to HD}}\sum _j{w}_{ij}^2{r}_j^{COMB}(t-\mathrm{\Delta }t)$(1)

where $h_i^{HD}(t)$ is the activation level of head direction cell i at time t. The term $-h_i^{HD}(t)$ represents an exponential decay term, whereby in the absence of presynaptic input, the activation of head direction cell i will fall to zero as determined by the time constant ${\tau }^{HD}$. The term $e_i(t)$ represents external visual input to head direction cell $i$ at the time $t$. This input term is used to generate an initial packet within the head direction cell layer. It is determined by a Gaussian profile described in pre-wired model simulation protocol (section 1.3). The term $-\frac{1}{N^{HD}}\sum _j{\tilde{w}}^{HD}{r}_j^{HD}(t)$ represents inhibitory feedback within the head direction cell layer, summed over all presynaptic head direction cells j. ${\tilde{w}}^{HD}$ is a global constant describing the effect of inhibitory interneurons, and $N^{HD}$ is the total number of head direction cells in the layer. The term $\frac{{\varphi }_2}{C^{COMB\to HD}}\sum _j{w}_{ij}^2{r}_j^{COMB}(t-\mathrm{\Delta }t)$ represents excitatory input from the COMB cell layer back onto the HD layer, summed across all presynaptic COMB cells connected to cell i. It consists of $r_j^{COMB}(t-\mathrm{\Delta }t)$, the firing of presynaptic COMB cell j at time $t-\mathrm{\Delta }t$, where $\mathrm{\Delta }t$ represents axonal conduction delay, and $w_{ij}^2$, the synaptic weight from presynaptic COMB cell j to postsynaptic HD cell i. It is important to note that neither the $w^2$ weights, nor any other weights, are updated during simulation, and thus remain at the same initial values throughout. The entire term is scaled by the factor $\frac{{\varphi }_2}{C^{COMB\to HD}}$, which controls the overall strength of COMB synapses onto the head direction cell layer where ${\varphi }_2$ is a constant and $C^{COMB\to HD}$ is the number of synapses each postsynaptic head direction cell receives from presynaptic COMB cells.

The firing rate r${ }_i$${ }^{\mathrm{H}\mathrm{D}}$(t) of head direction cell i at time t is a sigmoid function of the activation level of cell i

(2) $r_i^{HD}(t)=\frac{1}{1+e^{-2{\beta }^{HD}(h_i^{HD}(t)-{\alpha }^{HD})}}$(2)

where ${\alpha }^{HD}$ and ${\beta }^{HD}$ are the sigmoid threshold and slope, respectively. As this model is firing-rate based, no individual action potentials are simulated, and neuronal activity is represented as a spontaneous average firing rate.

The equations governing the activation level h${ }_i$${ }^{COMB}$ of ROT-COMB cells and NOROT-COMB cells indexed by i differ slightly. The ROT-COMB cell activation is given by

(3) $\begin{array}{rl}& {\tau }^{COMB}\frac{d{h}_i^{COMB}(t)}{dt}=-h_i^{COMB}(t)\\ & -\frac{1}{N^{COMB}}\sum _j{\tilde{w}}^{COMB}{r}_j^{COMB}(t)\\ & +\frac{{\varphi }_1}{C^{HD\to COMB}}\sum _j{w}_{ij}^1{r}_j^{HD}(t-\mathrm{\Delta }t)\\ & +\frac{{\varphi }_3}{C^{ROT\to COMB}}\sum _j{w}_{ij}^3{r}_j^{ROT}(t)\end{array}$(3)

whilst NOROT-COMB activation is given by

(4) $\begin{array}{rl}& {\tau }^{COMB}\frac{d{h}_i^{COMB}(t)}{dt}=-h_i^{COMB}(t)\\ & -\frac{1}{N^{COMB}}\sum _j{\tilde{w}}^{COMB}{r}_j^{COMB}(t)\\ & +\frac{{\varphi }_1}{C^{HD\to COMB}}\sum _j{w}_{ij}^1{r}_j^{HD}(t-\mathrm{\Delta }t)\\ & +\frac{{\varphi }_4}{C^{NOROT\to COMB}}\sum _j{w}_{ij}^4{r}_j^{NOROT}(t)\end{array}$(4)

where $h_i^{COMB}(t)$ is the activation level of COMB cell $i$ at time $t$. The term $-h_i^{COMB}(t)$ represents exponential decay, where in the absence of any input the activation of COMB cell $i$ will fall to zero as determined by the time constant ${\tau }^{COMB}$. The term $-\frac{1}{N^{COMB}}\sum _j{\tilde{w}}^{COMB}{r}_j^{COMB}(t)$ represents inhibitory feedback within the COMB cell layer, summed over all presynaptic COMB cells j. ${\tilde{w}}^{COMB}$ is a global constant describing the effect of inhibitory interneurons, and $N^{COMB}$ is the total number of COMB cells in the layer. The term $\frac{{\varphi }_1}{C^{HD\to COMB}}\sum _j{w}_{ij}^1{r}_j^{HD}(t-\mathrm{\Delta }t)$ represents excitatory input from the HD cell layer back onto the COMB layer, summed across all presynaptic HD cells connected to cell i. It consists of $r_j^{HD}(t-\mathrm{\Delta }t)$, the firing of presynaptic HD cell j at time $t-\mathrm{\Delta }t$, where $\mathrm{\Delta }t$ represents an axonal conduction delay, and $w_{ij}^1$, the synaptic weight from presynaptic HD cell j to postsynaptic COMB cell i. It is scaled by the factor $\frac{{\varphi }_1}{C^{HD\to COMB}}$, which controls the overall strength of HD synapses onto the COMB cell layer, where ${\varphi }_1$ is a constant and $C^{HD\to COMB}$ is the number of synapses each postsynaptic COMB cell receives from presynaptic HD cells.

The final term in the activation equation differs for NOROT-COMB and ROT-COMB cells. For ROT-COMB cells it is $\frac{{\varphi }_3}{C^{ROT\to COMB}}\sum _j{w}_{ij}^3{r}_j^{ROT}(t)$, where $\sum _j{w}_{ij}^3{r}_j^{ROT}(t)$ is the summed excitatory input from presynaptic ROT cells indexed by $j$ at time $t$, and $w_{ij}^3$ is the synaptic weight from presynaptic ROT cell j to postsynaptic COMB cell i. This is scaled by the factor $\frac{{\varphi }_3}{C^{ROT\to COMB}}$ where ${\varphi }_3$ is a constant and $C^{ROT\to COMB}$ is the number of synapses each postsynaptic COMB cell receives from presynaptic ROT cells. These terms are similar for NOROT-COMB cells, with the substitution of NOROT for ROT and ${\varphi }_4$ for ${\varphi }_3$. In these simulations, the number of ROT cells $N^{ROT}=1$ and the number of NOROT cells $N^{NOROT}=1$, and the synaptic weights $w_{ij}^3$ and $w_{ij}^4$ are both set to 1. As such the terms $\frac{{\varphi }_3}{C^{ROT\to COMB}}\sum _j{w}_{ij}^3{r}_j^{ROT}(t)$ and $\frac{{\varphi }_4}{C^{NOROT\to COMB}}\sum _j{w}_{ij}^4{r}_j^{NOROT}(t)$ can be simplified to ${\varphi }_3{r}^{ROT}(t)$ and ${\varphi }_4{r}^{NOROT}(t)$ respectively.

Note, then, that the activations of ROT-COMB and NOROT-COMB cells differ only by whether they receive input from ROT or NOROT cells. In all other regards, i.e., where inhibition is concerned, COMB cells function as a homogenous population.

Similar to HD cell firing rates, the firing rate r${ }_i$${ }^{\mathrm{C}\mathrm{O}\mathrm{M}\mathrm{B}}$(t) of COMB cell i at time t is a sigmoid function of the activation level of cell i

(5) $r_i^{COMB}(t)=\frac{1}{1+e^{-2{\beta }^{COMB}(h_i^{COMB}(t)-{\alpha }^{COMB})}}$(5)

Note that both the sigmoid thresholds ${\alpha }^{COMB}$ and ${\alpha }^{HD}$, as well as slopes ${\beta }^{COMB}$ and ${\beta }^{HD}$, may differ in value between COMB and HD cell layers.

In this pre-wired network, all the connectivity is set before simulation starts and does not change over the course of simulation. $w^1$ connectivity from HD to COMB cells is set-up in two distinct stages. Firstly, non-offset connectivity is set-up from HD to NOROT COMB cells. This is calculated as a Gaussian configuration, with the weight $w_{ij}^1$ for the synapse onto postsynaptic NOROT-COMB cell $i$ from presynaptic head direction cell $j$ being

(6) $w_{ij}^1=e^{-{(s_{ij}^{HDCOMB})}^2/2{({\sigma }^{HDCOMB})}^2}$(6)

where ${\sigma }^{HDCOMB}$ is the standard deviation of the Gaussian profile. $s_{ij}^{HDCOMB}$ is the difference between the preferred head directions $s_i$ and $s_j$ of the post- and pre-synaptic cells $i$ and $j$ respectively. It is given by

(7) $s_{ij}^{HDCOMB}=MIN(\left|x_i-x_j\left|,360-\right|x_i-x_j\right|)$(7)

where $x_i$ and $x_j$ are the preferred head directions of cells $i$ and $j$ respectively. This creates a wrap-around effect, with weights between head direction cells as a population remaining continuous across the $360/0{}^{\circ }$divide.

Offset connectivity from HD to ROT-COMB cells is calculated similarly, with the weight $w_{ij}$ for the synapse onto postsynaptic ROT-COMB cell $i$ from presynaptic head direction cell $j$ being calculated as in equation 6. However, the difference between preferred head directions of post- and pre-synaptic cells $i$ and $j$ now includes an additional offset $O$ and is calculated as

(8) $s_{ij}^{HDCOMB}=MIN(\left|x_i-(x_j+O)\left|,360-\right|x_i-(x_j+O)\right|)$(8)

where the offset $O$ is calculated, based on the target speed $V{}^{\circ }/s$ and the axonal conduction delay $\mathrm{\Delta }t$, as

(9) $O=V\mathrm{\Delta }t$(9)

This effectively causes presynaptic HD cells to be connected to an offset region of the ROT-COMB cell layer.

$w^2$ connectivity from COMB to HD cells is pre-wired in an identical fashion to $w^1$ connectivity. The sole difference is that the perspective is now switched, with the HD layer becoming post-synaptic cells indexed by $i$ and COMB cells are presynaptic cells $j$. The resulting combination of $w^1$ and $w^2$ connectivity means that individual HD cells both project symmetrically back onto themselves via NOROT-COMB cells, as well as asymmetrically back onto a region of the HD layer with an offset of $2O$ via ROT-COMB cells.

As this model makes use of a single ROT and a single NOROT cell, both COMB cell subtypes receive input from their resective AHV cell weighted directly by the parameters ${\varphi }_3$ for the ROT cell onto ROT-COMB cells and ${\varphi }_4$ for the NOROT cell onto NOROT-COMB cells.

The parameter values used for the pre-wired model are as shown in Table 1 unless stated otherwise.

Table 1. Network parameter values used during simulation of the pre-wired model. All pre-wired simulations, except where explicitly stated, are run using these parameters. All measures of time are given in seconds.

Pre-wired Network Parameters
No. Head Direction Cells	500
No. ROT-COMB Cells	500
No. NOROT-COMB cells	500
Simulation Time	$4.1s$
$C^{COMB\to HD}$	1000
$C^{HD\to COMB}$	500
ϕ1	700.0
ϕ2	4500.0
${\sigma }^{HDCOMB}$	${20.0}^{\circ }$
${\sigma }^{COMBHD}$	${20.0}^{\circ }$
No. ROT Cells	1
No. NOROT Cells	1
$C^{ROT\to COMB}$	1
$C^{NOROT\to COMB}$	1
ϕ3	80.0
ϕ4	80.0
${\alpha }^{HD}$	0.0
${\beta }^{HD}$	0.2
${\alpha }^{COMB}$	16.0
${\beta }^{COMB}$	0.3
${\tau }^{HD}$	$0.0001s$
${\tau }^{COMB}$	$0.0001s$
${\tilde{w}}^{HD}$	0.2
${\tilde{w}}^{COMB}$	0.35
$\delta t$	$0.00001s$
${\lambda }^{HD}$	2.0
${\sigma }^{HD}$	${20.0}^{\circ }$
Target Velocity V	$180{}^{\circ }/s$
Axonal delay $\mathrm{\Delta }t$	$0.01s$

The differential equations given for this model cannot be solved analytically. They are solved numerically by a Forward Euler finite difference scheme with a timestep size of $\delta t$ (as specified in Table 1).

Pre-wired model simulation protocol

At the beginning of the simulation, the firing rates and activations of all cells are set to zero. The connectivity from HD to COMB, from COMB to HD, and from AHV cells to COMB is set. An external visual input, represented as the term $e_i$, is applied to each cell $i$ in the HD layer for $100ms$ to initialize network activity. $e_i$ is calculated for a given postsynaptic HD cell $i$ according to the Gaussian function

(10) $e_i={\lambda }^{HD}{e}^{-{(s_i^{HD})}^2/2{({\sigma }^{HD})}^2}$(10)

where ${\lambda }^{HD}$ is a scaling factor determining the strength of this input to head direction cell $i$ and ${\sigma }^{HD}$ is the standard deviation of the Gaussian profile. $s_i^{HD}$ is the difference between the true head direction $x$ and the preferred direction $x_i$ of postsynaptic head direction cell $i$. It is given by

(11) $s_i^{HD}=MIN(\left|x_i-x\left|,360-\right|x_i-x\right|)$(11)

which creates a wrap-around effect, with response profiles of head direction cells as a population remaining continuous across the $360/0{}^{\circ }$divide. $e_i$ thus generates in the HD layer a Gaussian activity packet centred on location $x$. During this period, HD activations $h_i^{HD}$ and firing rates $r_i^{HD}$ are updated as in equations 1 and 2. COMB activations $h_i^{COMB}$ and firing rates $r_i^{COMB}$ are also updated as in equations 3, 4, and 5, with $r^{NOROT}=1$ and $r^{ROT}=0$.

After this initialization period, $e_i$ is removed. There then follows a $1s$ period of zero head rotation. All activations and firing rates continue to update, and AHV cells remain set with $r^{NOROT}=1$ and $r^{ROT}=0$.

Next, there is a $2s$ period of head rotation at $V{}^{\circ }/s$. AHV cell firings are changed to reflect the fact that head rotation is occurring, with $r^{NOROT}=0$ and $r^{ROT}=1$. HD and COMB cell activations and firing rates continue to update.

Finally, there is another $1s$ period of zero head rotation, in which HD and COMB cell activations and firings continue to update and AHV cell firings are set to reflect zero head rotation with $r^{NOROT}=1$ and $r^{ROT}=0$.

Self-organizing model description

Given that pre-wired simulations demonstrate the core functionality of the hypothesized system with idealized connectivity, a next step is to check whether such a system can self-organize through learning. Can the self-organized model, in which the synaptic connections have been set-up by learning, perform accurate path integration?

In the self-organizing model, we do not pre-categorize COMB cells into ROT-COMB and NOROT-COMB types. Instead, we initially lump all COMB cells into a single undifferentiated category, denote both types of AHV cells as ROT cells, and then connect all ROT cells to all COMB cells. The $w^3$ connections from ROT cells to COMB cells are then set-up during learning, which will determine the ultimate response selectivity of individual COMB cells to ROT cells which are active when the head is rotating and ROT cells which are active when the head is not rotating. Can ROT-COMB and NOROT-COMB cells be formed without this distinction being imposed on the model by artificial pre-wiring of the synaptic connections?

The self-organizing variant is very similar to the pre-wired model. The principal difference is that connectivity is not pre-wired. All connections (HD-COMB, COMB-HD, ROT-COMB) are initialized with random values. During simulation, all weight profiles self-organize via Hebbian learning rules over the course of model training in the light, detailed below. The simulation protocol is then identical to the pre-wired model.

Crucially, HD to COMB connectivity is diluted. That is, each COMB cell receives projections from a limited number of HD cells. In the simulations detailed below, each of the $1000$ COMB cells receives from $25$ of the $500$ HD cells, giving $5\mathrm{\%}$ connectivity. This is in order to avoid the effects of continuous transformation (CT) learning (Stringer et al. 2006). In brief, the CT effect is a result of the fact that the HD activity packet moves with spatiotemporal continuity during training. Say that an arbitrary pattern of COMB cell firing is driven up by an HD activity packet at a given position. Connections will be strengthened, by Hebbian learning, between all active HD cells and all active COMB cells. If at the next timestep, the next HD activity packet has enough overlapped with that of the previous timestep, the same or very similar pattern of COMB cells will be driven up by virtue of the fact that a similar set of HD cells is firing. However, some additional HD cells will become active, and their connections onto active COMB cells will be strengthened by Hebbian learning. Over multiple iterations of this process, it is possible for connections to be strengthened between a particular subset of COMB cells and all HD cells, thus eliminating COMB specificity to HD. By reducing the number of HD cells from which a COMB cell receives connections, the CT effect is avoided.

Also, present in the self-organizing model is feed-forward inhibition from the visual input to the HD layer during training. This is to ensure that the effect of excitatory visual input during training does not flood the network with activation when combined with HD-COMB reciprocal connectivity.

Unlike the pre-wired model, no a priori distinction is made between ROT-COMB and NOROT-COMB cells. Instead, the COMB layer receives connections from a layer of 1000 ROT cells, half of which are active during head rotation, and half of which are active when the head is stationary.

Self-organizing model equations

The HD layer activation equation now incorporates feed-forward inhibition from visual input. This maintains a similar level of network activation across both training in the light and testing in the dark and is represented by the constant $I^{FF}$. The activation level h${ }_i$${ }^{HD}$ of head direction cell i is therefore now given by

(12) $\begin{array}{rl}& {\tau }^{HD}\frac{d{h}_i^{HD}(t)}{dt}=-h_i^{HD}(t)+e_i(t)\\ & -\frac{1}{N^{HD}}\sum _j{\tilde{w}}^{HD}{r}_j^{HD}(t)\\ & +\frac{{\varphi }_2}{C^{COMB\to HD}}\sum _j{w}_{ij}^2(t)r_j^{COMB}(t-\mathrm{\Delta }t)\\ & -I^{FF}\end{array}$(12)

where $-I^{FF}$ represents feedforward inhibition from the visual input.

The firing rate r${ }_i$${ }^{\mathrm{H}\mathrm{D}}$(t) of head direction cell i at time t is a sigmoid function of the activation level of cell i given by equation 2.

All COMB cells receive from all ROT cells, via $w^3$ weights. The activation level h${ }_i$${ }^{COMB}$ of COMB cell $i$ is given by

(13) $\begin{array}{rl}& {\tau }^{COMB}\frac{d{h}_i^{COMB}(t)}{dt}=-h_i^{COMB}(t)\\ & -\frac{1}{N^{COMB}}\sum _j{\tilde{w}}^{COMB}{r}_j^{COMB}(t)\\ & +\frac{{\varphi }_1}{C^{HD\to COMB}}\sum _j{w}_{ij}^1(t)r_j^{HD}(t-\mathrm{\Delta }t)\\ & +\frac{{\varphi }_3}{C^{ROT\to COMB}}\sum _j{w}_{ij}^3(t)r_j^{ROT}(t)\end{array}$(13)

The final term is $\frac{{\varphi }_3}{C^{ROT\to COMB}}\sum _j{w}_{ij}^3(t)r_j^{ROT}(t)$, which represents ROT cell influence on the COMB cell layer, where $\sum _j{w}_{ij}^3(t)r_j^{ROT}(t)$ is the summed firing of presynaptic ROT cells indexed by $j$ at time $t$ multiplied by the weight $w_{ij}^3(t)$ from presynaptic ROT cell $j$ to postsynaptic COMB cell $i$ at time t. This is scaled by the factor $\frac{{\varphi }_3}{C^{ROT\to COMB}}$ where ${\varphi }_3$ is a constant, and $C^{ROT\to COMB}$ is the number of synapses each postsynaptic COMB cell receives from presynaptic ROT cells.

Similarly, to HD cells, the firing rate r${ }_i$${ }^{\mathrm{C}\mathrm{O}\mathrm{M}\mathrm{B}}$(t) of COMB cell i at time t is a sigmoid function of the activation level of cell i given by equation 5

Synaptic weights are updated during training based on Hebbian learning rules. These rules take the same form regardless of whether connections are HD-COMB, COMB-HD, or ROT-COMB. The weight between presynaptic HD cell $j$ and postsynaptic COMB cell $i$ is updated according to a local associative Hebbian learning rule

(14) $\frac{d{w}_{ij}^1(t)}{dt}=k{r}_i^{COMB}(t)r_j^{HD}(t-\mathrm{\Delta }t)$(14)

where $w_{ij}^1$ is the synaptic weight from presynaptic HD cell $j$ with firing rate $r_j^{HD}$ onto postsynaptic COMB cell $i$ with firing rate $r_i^{COMB}$, and $k$ is the learning rate constant, which determines the speed of weight change.

Similarly, the weight between presynaptic COMB cell $j$ and postsynaptic HD cell $i$ is updated by the equation

(15) $\frac{d{w}_{ij}^2(t)}{dt}=k{r}_i^{HD}(t)r_j^{COMB}(t-\mathrm{\Delta }t)$(15)

Finally, the weight between presynaptic ROT cell $j$ and postsynaptic COMB cell $i$ is updated by the equation

(16) $\frac{d{w}_{ij}^3(t)}{dt}=k{r}_i^{COMB}(t)r_j^{ROT}(t)$(16)

At each timestep all synaptic weight vectors are also normalized after updating. This is achieved by ensuring that the square root of the sum of the squares of afferent weights, for a particular postsynaptic cell i, is limited to 1. The weight normalization equation for HD to COMB weights is given as

(17) $\sqrt{\sum _j{(w_{ij}^1(t))}^2}=1$(17)

and for COMB to HD weights

(18) $\sqrt{\sum _j{(w_{ij}^2(t))}^2}=1$(18)

and for ROT to COMB weights

(19) $\sqrt{\sum _j{(w_{ij}^3(t))}^2}=1$(19)

The parameter values used for the self-organizing model are given in Table 2 unless stated otherwise.

Table 2. Effect of varying values of $\mathrm{\Delta }t$ for constant ${\tau }^{HD/COMB}$ on an average time interval between HD activity packet shifts, between COMB shifts, and between HD and COMB shifts. All measures of time are given in seconds.

Effect of varying $\mathrm{\Delta }t$ for constant $\tau$ on system dynamics
${\tau }^{HD/COMB}$	$\mathrm{\Delta }t$	HD-HD Interval	COMB-COMB Interval	HD-COMB Interval
0.0001	0.0025	0.005	0.005	0.0025
0.0001	0.005	0.01	0.01	0.005
0.0001	0.0075	0.015	0.015	0.0075
0.0001	0.01	0.02	0.0201	0.0099
0.0001	0.025	0.05	0.05	0.025
0.0001	0.05	0.1	0.1	0.05

Self-organizing model simulation protocol

At the beginning of the simulation, the firing rates and activations of all cells are set to zero. Connectivity from HD to COMB, from COMB to HD and from ROT to COMB is set. In these self-organizing simulations, all weights are initialized with random values. This is achieved using the GNU C Library’s srand48() and drand48(). Together, these functions generate non-negative random numbers drawn from a uniform distribution in the interval $[0.0,1.0]$. Following random initialization of weight values, all synaptic weight vectors are normalized. In the case of diluted HD to COMB connectivity, $25$ connections for each COMB cell are randomly selected from the $500$ possible presynaptic HD cells using the gsl_ran_choose() function.

Following this initialization, training begins. A single training epoch is comprized of rotating and stationary periods. During the rotation period, the network undergoes a period of an HD activity packet being imposed by external visual input (as described for the pre-wired simulations) and shifted through the HD ring at the target velocity $V{}^{\circ }/s$. The firing rates of half of the ROT cells are set to 1, whilst the other half of the ROT cells have firing rates of 0. Following the rotation period, ROT cells switch states and there follows a stationary period. During this latter period, the activity packet moves to be centred on each HD cell in turn for a fixed period of time. The amount of time spent during the rotation period and the stationary period is kept equal, such that the amount of time spent rotating is 500 (the number of discrete HD cell locations) times the amount of time spent stationary at each location. This is to ensure equal training with the agent rotating and stationary. This process constitutes a single epoch of training and can be repeated as much as is required.

After training, the 4-second testing period begins. External visual input and feed-forward inhibition are both removed. In the first second of testing, the HD packet remains where it ended up after training, and those ROT cells active during stationary training are active. This holds the HD activity packet at the current location. During the next 2 s, ROT cells switch states (i.e., those signals rotation become active, whilst stationary ROT cells are quiescent), and the motion of the HD activity packet is observed. In the final second of testing, ROT cells switch states once more, and the HD activity packet is held at its final location.

General operating principles

A diagram of the pre-wired and self-organizing models’ operation is given in Figure 2. When the agent’s head is not rotating, NOROT-COMB cells (which are generated either by pre-wiring or via self-organization), responding to a combination of zero angular head velocity and the current head direction, will become active. These cells project directly back to the HD cell from which they receive. This non-offset reciprocal connectivity across two layers will allow the HD activity packet to be sustained in its current location. Crucially, this connectivity will not remain active during head rotation, as the angular head velocity cells will no longer signal a lack of head rotation, and thus NOROT-COMB cells will become quiescent.

Figure 2. Operation of the pre-wired and self-organizing models/HD activity is both sustained and shifted by reciprocal connectivity with a layer of combination (COMB) cells. This connectivity contains a conduction delay $\mathrm{\Delta }t$. COMB cells respond to a particular combination of HD and angular head velocity (AHV). During no head rotation (a), HD activity is stabilized through reciprocal connectivity with COMB cells representing the combination of current head direction and no angular head velocity (AHV = 0). When the head is rotating (b), activity is shifted through the HD layer via connectivity with different COMB cells representing the conjunction of current head direction with angular head velocity V (AHV = V, where V is some non-zero value). Connectivity from these COMB cells back onto the HD layer projects with an offset from the original HD activity of $2\mathrm{\Delta }tV$; causing HD cell activity to accurately track true head direction given current angular head velocity V.

When the agent rotates its head at arbitrary non-zero velocity $V{}^{\circ }/s$, ROT-COMB cells will become active. These cells respond to a combination of current head direction and angular head velocity $V{}^{\circ }/s$. They project back onto the HD layer with an offset based on the target velocity and the conduction delay $\mathrm{\Delta }t$. This offset connectivity means that a given HD cell will project, via the ROT-COMB cells, most strongly to the HD cell which ought to be active $2\mathrm{\Delta }t$ later, given the total conduction delay $2\mathrm{\Delta }t$ and the target velocity $V{}^{\circ }/s$. It should be noted that the effective total conduction delay, in this case, will be $2\mathrm{\Delta }t$, as the HD layer projects back onto itself via two sets of connections (HD to COMB and COMB to HD). This will act to shift the activity through the HD layer to the correct location, allowing the system to accurately track true head direction given current angular head velocity. Effectively, this HD to ROT-COMB to HD connectivity allows the system to know where a final HD activity packet ought to be, given a particular starting location, following head rotation at a specific angular head velocity over a known time interval.

Data analysis

In order to calculate the speed of motion of the HD layer activity packet, the centre of mass of HD layer firing rates (i.e., the location of HD activity packet, also known as the population or Pvector) is computed at each forward Euler timestep. This was calculated as a circular mean, using an established population vector scheme (Georgopoulos et al. 1986; Song and Wang 2005) as follows

(20) ${\theta }_{pop}(t)=arctan\left(\frac{\sum _i{r}_{i}sin({\theta }_i)}{\sum _i{r}_{i}cos({\theta }_i)}\right)$(20)

where ${\theta }_i$ is the preferred direction of HD cell $i$ with firing rate $r_i$.

To give a value for the population vector direction in the range [$0-360{}^{\circ }$], the following corrected formula is used instead, with this functionality being captured by the MATLAB function atan2

$arctan\left(\frac{\sum _i{r}_{i}sin({\theta }_i)}{\sum _i{r}_{i}cos({\theta }_i)}\right)\mathrm{i}\mathrm{f}\sum _i{r}_{i}sin({\theta }_i)>0,\sum _i{r}_{i}cos({\theta }_i)>0;$

${\theta }_{pop}(t)=arctan\left(\frac{\sum _i{r}_{i}sin({\theta }_i)}{\sum _i{r}_{i}cos({\theta }_i)}\right)+180\mathrm{i}\mathrm{f}\sum _i{r}_{i}cos({\theta }_i)<0;$

(21) $arctan\left(\frac{\sum _i{r}_{i}sin({\theta }_i)}{\sum _i{r}_{i}cos({\theta }_i)}\right)+360\mathrm{i}\mathrm{f}\sum _i{r}_{i}sin({\theta }_i)<0,\sum _i{r}_{i}cos({\theta }_i)>0;$(21)

For the pre-wired model, the COMB cells are split into ROT-COMB and NOROT-COMB populations, with each population each having its own distribution of cell preferred firing directions in the interval [$0-360{}^{\circ }$]. The above Pvector calculation is therefore applied separately to ROT-COMB and NOROT-COMB cells. However, for ease of visualization, Pvectors are collated and thus will show NOROT-COMB Pvectors during stationary phases, and ROT-COMB Pvectors during rotation. Pvectors are also used to calculate path integration speed, which was calculated by noting the total distance travelled during the rotation phase of simulation and dividing by rotation time.

Results

Pre-wired model is capable of extremely accurate path integration

Figure 3 summarizes HD layer packet speed during simulated path integration with a variety of target velocities in the range [$0-360{}^{\circ }/s$]. Parameter values are given in Table 1. We see clearly that the path integration velocity attainable by the model is extremely close to the target velocity. Indeed, this value is in the region of $99.5\mathrm{\%}$ for all target velocities in the range tested.

Figure 3. Simulation results of the pre-wired network model. Path integration speed is shown as a percentage of the target for a range of target velocities, with parameters as given in Table 1. Model performance is excellent (observed speed $>99\mathrm{\%}$ of target) at all target velocities tested.

Effect of time constant

A key factor in the accuracy of path integration speed (see Page et al., 2015), is the neuronal time constant ${\tau }^{HD}$ (and in this model ${\tau }^{COMB}$). For this reason, an investigation was made into the effects of time constant across values of ${\tau }^{HD}$ and ${\tau }^{COMB}$ in the range [$0.0001s$,$0.1s$]. The results of this investigation are given in Figure 4. In all cases, target velocity was $V=180{}^{\circ }/s$, conduction delay was fixed at $\mathrm{\Delta }t=0.01$, and ${\tau }^{HD}={\tau }^{COMB}$. As expected, we see that increasing rise time via increasing the value of both ${\tau }^{HD}$ and ${\tau }^{COMB}$ reduces packet speed during path integration.

Figure 4. Effect of neuronal time constant on path integration accuracy.

Speeds during path integration are shown for a range of values of ${\tau }^{HD}={\tau }^{COMB}$ with a target velocity $V=180{}^{\circ }/s$. As ${\tau }^{HD}$ and ${\tau }^{COMB}$ increase, so too does the effect of rise time, and thus path integration speed is reduced.

Effect of conduction delay

Another key factor in path integration speed is the conduction delay $\mathrm{\Delta }t$. A variety of different conduction delays were used in this model. Figure 5 shows the summary results from simulations run with the parameters given in Table 1 and $\mathrm{\Delta }t$ varying in the range [$0.001s$,$0.05s$]. In all simulations ${\tau }^{HD}={\tau }^{COMB}=0.001s$ and $V=180{}^{\circ }/s$. We see that shorter conduction delays are associated with reduced packet speed, whilst longer conduction delays show packet speeds closer to the target. This is because a constant error introduced by neuronal rise time will be proportionally more severe in the context of a smaller conduction delay rather than a larger one (Page et al. 2015).

Figure 5. Effect of conduction delay on path integration accuracy.

Speeds during path integration are shown for a range of values of conduction delay $\mathrm{\Delta }t$ with a target velocity of $V=180{}^{\circ }/s$ and a fixed neuronal time constant of ${\tau }^{HD}={\tau }^{COMB}=0.001s$. A constant error introduced by rise time will be proportionally more severe in the context of shorter conduction delays. Thus, the slowing effect on path integration of rise time is increased for small $\mathrm{\Delta }t$ and decreased for large $\mathrm{\Delta }t$.

The results so far demonstrate model functionality when all HD to COMB $w^1$ and COMB to HD $w^2$ connections involve the same axonal conduction delay $\mathrm{\Delta }t$. However, it is highly unlikely that connectivity within the HD system or indeed the brain as a whole, is associated with delays of a single uniform value. Simulations in which all $w^1$ and $w^2$ synapses have their own value of $\mathrm{\Delta }t$, drawn randomly from a uniform distribution in the range [$0.0001s$,$0.1s$], were therefore carried out. In these simulations, the effect of time constant was found to be similar to that with a single value of $\mathrm{\Delta }t$, and these results are pictured for a target velocity of $V=180{}^{\circ }/s$ in Figure 6. One crucial difference, evident in this Figure, is that path integration was consistently faster with a range of conduction delay values than with a single value. Although the speed when ${\tau }^{HD}={\tau }^{COMB}=0.0001s$ was similar ($179.43{}^{\circ }/s$ in these simulations, compared with ${179.73}^{\circ }/s$ in simulations with a single $\mathrm{\Delta }t$ value), comparison of Figures 4 and 6 reveals that large time constants are less detrimental to path integration accuracy in a model with a distribution of $\mathrm{\Delta }t$ values. This suggests that varying values of conduction delay are a partial solution to the issue of reduced path integration accuracy with larger neuronal time constants, possibly helping to ameliorate the effect of rise time.

Figure 6. Effect of neuronal time constant on a model with a distribution of conduction delay values in the range [$0.0001s$, $0.1s$].

The model utilizing a distribution of conduction delays maintains significantly better path integration accuracy across the range of $\tau$ values when compared to a model with a single value of $\mathrm{\Delta }t$. Identically to Figure 4, we see that increased rise time via increasing ${\tau }^{HD}={\tau }^{COMB}$ reduces path integration accuracy.

Periodic behaviour

In simulations with a single value of $\mathrm{\Delta }t$, examining the motion of HD and COMB activity packets during path integration more closely reveals that both HD and COMB layer packet motion is not smooth, but consists of periodic movements. This is likely an effect of having two linked layers separated by a single value of conduction delay: an HD activity packet will shift itself, via the COMB layer, after $2\mathrm{\Delta }t$. This intuition is confirmed by Figure 7, which shows packet location from 1 s to 1.4 s of simulation with both a single fixed value of $\mathrm{\Delta }t=0.005s$ (left) and a distribution of $\mathrm{\Delta }t$ values in the range [$0.001s$,$0.01s$]. With a single value, packet motion is not continuous, but instead shows a periodic pattern of shifting. In contrast, a distribution of conduction delays leads to both HD (blue) and COMB (red) activity packets moving continuously without periodic jumps.

Figure 7. Periodic motion with a single $\mathrm{\Delta }t$ value.

Position (PVector) of both HD (blue) and COMB (red) activity packets from seconds 1 to 1.4 of the simulation. Periodic motion is observed for a single value of $\mathrm{\Delta }t=0.005s$ (left) but is not observed when using a distribution of $\mathrm{\Delta }t$ values in the range [$0.001s$,$0.01s$] (right).

Manipulating $\mathrm{\Delta }t$ and ${\tau }^{HD/COMB}$ provides more insight into the specific nature of this periodic phenomenon, and verifies that the pre-wired model is operating as hypothesized. Tables 2 and 3 show, respectively, the effects of varying $\mathrm{\Delta }t$ for constant ${\tau }^{HD/COMB}$ and varying ${\tau }^{HD/COMB}$ for constant $\mathrm{\Delta }t$. Periodic behaviour matches the operating principles of the model: HD shifts are separated from one another by $2\mathrm{\Delta }t$, which is also the separation between COMB shifts. HD shifts are separated from COMB shifts by $\mathrm{\Delta }t$. This confirms that the model operates as hypothesized.

Table 3. Effect of varying values of ${\tau }^{HD/COMB}$ for constant $\mathrm{\Delta }t$ on an average time interval between HD activity packet shifts, between COMB shifts, and between HD and COMB shifts. All measures of time are given in seconds.

Effect of varying $\tau$ for constant $\mathrm{\Delta }t$ on system dynamics
${\tau }^{HD/COMB}$	$\mathrm{\Delta }t$	HD-HD Interval	COMB-COMB Interval	HD-COMB Interval
0.0001	0.05	0.1	0.1	0.05
0.001	0.05	0.0999	0.0999	0.05002
0.0025	0.05	0.0999	0.0999	0.05028
0.005	0.05	0.1047	0.1047	0.05053
0.0075	0.05	0.1088	0.1091	0.053
0.01	0.05	0.1136	0.1146	0.0549

Varying $\mathrm{\Delta }t$ affects periodic HD and COMB layer packet shifts in a systematic manner. This was checked by re-running simulations with varying $\mathrm{\Delta }t$ values in the range [$0.001s$,$0.05s$], this time with $\tau$ reduced to a negligibly small value ${\tau }^{HD}={\tau }^{COMB}=0.0001s$ in order to see more clearly the effect of changing conduction delay without interference from neuronal rise time (Table 2). As would be expected given model operating principles, the inter-shift interval is always $2\mathrm{\Delta }t$ for HD-HD and COMB-COMB intervals and is always $\mathrm{\Delta }t$ for HD-COMB intervals.

Periodic behaviour was affected by rise time (Table 3), controlled by ${\tau }^{HD/COMB}$. Whilst periodic motion was still observed, the regular relationship between its timing and conduction delay broke down as rise time increased. This is due to increasing rise time causing neurons to start firing later relative to when they are first stimulated, with the intervals between HD-HD, COMB-COMB, and HD-COMB, shifts becoming larger than expected for a given value of $\mathrm{\Delta }t$. This clearly demonstrates rise time’s specific effect on path integration speed.

Self-organization

The self-organizing variant of the model was simulated with a target velocity of $180{}^{\circ }/s$ and all parameters as in Table 4. Figure 8 shows HD layer firing rates and position during testing. As we can see, the self-organizing network has made a good job of replicating the results of a pre-wired network, with approximately two stationary phases and a rotation phase, with the switch between the two phases governed by the activity of ROT cells. This demonstrates clearly that model functionality can be self-organized, with a single network being able to smoothly transition between holding and shifting an activity packet.

Table 4. Network parameter values used during self-organizing simulations. All self-organizing simulations, except where explicitly stated, are run using these parameters. All measures of time are given in seconds.

Self-organizing network parameters
No. Head Direction Cells	500
No. COMB Cells	1000
No. ROT cells	1000
Rotation Time (per training epoch)	$12.0s$
Stationary Time (per HD cell)	$0.024s$
No. Training Epochs	$50$
$C^{COMB\to HD}$	1000
$C^{HD\to COMB}$	25
ϕ1	2400.0
ϕ2	10000.0
$C^{ROT\to COMB}$	1000
ϕ3	5000.0
${\alpha }^{HD}$	0.0
${\beta }^{HD}$	0.2
${\alpha }^{COMB}$	113.0
${\beta }^{COMB}$	0.3
${\tau }^{HD}$	$0.0001s$
${\tau }^{COMB}$	$0.0001s$
${\tilde{w}}^{HD}$	1.0
${\tilde{w}}^{COMB}$	1000.0
$\delta t$	$0.00001s$
${\lambda }^{HD}$	450.0
${\sigma }^{HD}$	${20.0}^{\circ }$
$I^{FF}$	$420.0$
Target Velocity V	$180{}^{\circ }/s$
Learning rate k	$0.1$
Axonal delay $\mathrm{\Delta }t$	$0.01s$

Figure 8. Simulation results from the self-organizing model.

HD layer firing, with high firing indicated by black (top) and packet location (bottom) during the testing phase of simulation. Model performance is approximately correct, with a 2s period of packet motion during path integration, bookended by two 1 s periods of no packet motion whilst HD does not change.

The speed achieved during head rotation is $164.4{}^{\circ }/s$, which represents $91\mathrm{\%}$ of the target velocity. This is an improvement on previous self-organizing path integration models relying on within-layer recurrent collateral connections (Walters and Stringer 2010; Walters et al. 2013). This improvement in speed comes from the omission of within-layer recurrent collateral connectivity. However, the self-organizing model falls short of the ideal result achieved in the pre-wired model.

This difference must be related to differences in the weight structure that self-organizes. Inaccuracy will be introduced into the model if the wrong region of HD space is stimulated. This could well arise from poor COMB specificity to either HD or ROT or from poor specificity in the projections from COMB back onto HD.

To answer the question of specificity to head rotation, a ROT index statistic was constructed. For each COMB cell, the proportion of the total afferent ROT to COMB weights for that cell was calculated for the first 250 ROT cells (i.e., those ROT cells firing when the head is not rotating) and the second 250 ROT cells (i.e., those ROT cells firing when the head is rotating). This proportion is therefore in the range $[0,1]$. The ROT index is then calculated as the proportion for head rotation minus the proportion for no head rotation. Thus, each COMB cell has a ROT index ranging from $+1$ (all weight from ROT cells firing to head rotations) to $-1$ (all weight from ROT cells firing when the head is stationary). A ROT index of 0 would indicate no specificity to head rotation.

Figure 9 shows a histogram of ROT index values for the trained COMB cell layer. In total, 768 of the 1000 COMB cells ($76.8\mathrm{\%}$) had an absolute ROT index greater than 0.8. Of these cells, equal proportions were specific to head rotation ($37.5\mathrm{\%}$) or no head rotation ($37.5\mathrm{\%}$). However, the majority of the remaining $23.2\mathrm{\%}$ of COMB cells had a ROT index in the range $-0.6$ to $0.2$, indicating that COMB layer specificity to head rotation, whilst quite good, was not perfect.

Figure 9. Specificity of ROT to COMB synaptic weights.

In this histogram, a ROT index of $+1$ indicates that all afferent connection weight for that cell comes from ROT cells signalling head rotation, whilst a ROT index of $-1$ indicates that all afferent connection weight for that cell comes from ROT cells signalling a stationary head. ROT index values in between this range indicate mixed ROT cell specificity for a given COMB cell. This figure demonstrates that whilst the majority of COMB cells are specific to either head rotation or no head rotation, there is a noticeable proportion which are not.

To answer the question of COMB specificity to HD, Figure 10 shows a histogram of the number of COMB firing peaks, defined as contiguous periods of elevated COMB cell firing, as the model rotated through the full $360{}^{\circ }$ space of head directions during the test phase. Whilst many COMB cells either do not fire at all during the test phase or fire only once, we see a significant proportion of COMB cells show firing peaks twice or more during this same period. Given that the test phase consists of one complete rotation through the HD layer, COMB cells would be expected to fire just once during this period if specific to one HD.

Figure 10. Specificity of COMB cell response.

Histogram detailing the number of firing rate peaks (contiguous periods of elevated firing) for each cell in the COMB layer during testing as the model rotates through the full $360{}^{\circ }$ space of head directions. A notable proportion of cells has less specific responses to HD in that they display two or more firing peaks.

Conclusions

The major contribution of this paper is to offer a theoretical explanation for why there do not appear to be within-layer recurrent excitatory connections within distinct layers of neurons in the rat head direction system. Previous modelling studies have included such connections within the HD layer, to stabilize an activity packet in the absence of visual input (Walters and Stringer 2010; Walters et al. 2013). However, this always slows down path integration by around 20–30% and leads to an error in the representation of HD. The function of such a symmetric recurrent connectivity is merely to keep an activity packet rooted to the spot, as opposed to shifting the activity packet. The major problem that occurs is that the stabilizing recurrent connections cannot be switched off when the network needs to shift the packet at the correct speed during path integration of head rotation. The recurrent connections continue to exert their stabilizing influence, thus retarding the motion of the activity packet and introducing significant path integration error.

A more sophisticated learning rule, perhaps with an explicit axonal transmission delay, could self-organize asymmetrical recurrent connections tuned to a certain speed of head rotation. However, the core problem persists and can be stated more generally: any recurrent excitatory connections calibrated to a particular speed of head rotation will degrade path integration accuracy. Recurrent weight profiles can learn to respect at most one head velocity, and path integration will accquire inaccuracy during rotation at any other velocity.

This problem could be solved if there were multiple different sets of within-layer recurrent connections, each optimized for different angular head velocities, and each active only during appropriate kinds of head rotation. This is not possible within a biologically plausible network architecture. However, the novel contribution of this paper is to achieve something equivalent in a multi-layer form. Our model includes an additional layer of combination cells, which learn to respond to particular combinations of head direction and angular head velocity. Different subsets of combination cells become activated depending on the current angular head velocity of the simulated rat. Those cells responding to zero angular head velocity project back to the same HD layer location from which they receive input, whilst those cells responding to non-zero velocity project back to the HD layer with an offset.

By introducing an additional layer of COMB cells, with different subsets of COMB cells tuned to different angular head velocities, this effectively creates different subsets of bi-directional connections between the HD and COMB cell layers that may learn to implement path integration at corresponding velocities. In this manner, path integration input to the HD layer is always ”pure”, as the HD layer inputs signalling a single angular head velocity at any given time. For example, one subset of COMB cells and the corresponding subset of bi-directional connections between the HD and COMB cell layers may be tuned to no head rotation, and will be responsible for stabilizing the activity packet in the HD ring during periods of no head rotation in the dark. However, these COMB cells will respond specifically to zero angular head velocity, and will be inactive during any actual head rotation. Thus, the subset of stabilizing bi-directional connections will not interfere with path integration during head rotations.

This is why the simulations of path integration presented in this paper are much more accurate than our earlier models that relied on recurrent excitatory connections within the HD ring to stabilize the activity packet when the agent is stationary in the dark.

Discussion

This paper presents operational details and simulation results from pre-wired and self-organizing versions of a head direction cell path integration model which does not contain within-layer symmetric excitatory recurrent collateral weights. This work supports the assertion that a multi-layer system with bi-directional connections between layers exists principally for reasons of path integration speed (Page et al. 2015). The model in this paper consists of a ring of HD cells, reciprocally connected to a ring of COMB cells, with the latter also receiving input from ROT cells. All HD to COMB and COMB to HD connections contain an axonal conduction delay $\mathrm{\Delta }t$. COMB cells which receive input from ROT cells signalling no head rotation project back to the HD cell from which they maximally receive. These cells are responsible for sustaining a stationary HD layer activity packet when the head is not moving. COMB cells which receive input from ROT cells signalling head rotation at a target velocity project back onto the HD layer with an offset determined by the target velocity and the conduction delay. In this manner, direct HD to HD connections are no longer present, and instead, all HD cells influence one another via the COMB cell layer regardless of head motion.

Simulation results show that this model, when pre-wired with ideal connectivity, is capable of extremely accurate ($99\mathrm{\%}$) path integration, and this result is found to extend to a more biologically plausible situation in which conduction delays are randomly varied within a set range. Parameter exploration reveals that the operation of this model is much more accurate than previous self-organizing path integration models incorporating within-layer excitatory recurrent collateral connectivity (Walters and Stringer 2010; Walters et al. 2013). However, the effect of increasing rise time is, as with previous models, detrimental to path integration speed. This is because rise time causes neurons to fire slightly after when they ought to fire.

Interestingly, a periodic effect is seen in the model, whereby a fixed conduction delay value across all synapses causes activity in the model to shift in a stepwise manner, with the time for each step defined by conduction delay. Although it is shown that a distribution of conduction delay values, assumed to be more biologically plausible, eliminates this periodicity, the model does predict that a sufficiently narrow distribution of delays (i.e., one with low variance) in vivo would cause such periodicity.

It is also demonstrated that this kind of connectivity is capable of self-organizing, although path integration accuracy is slightly reduced in this case. What might be the reason for sub-optimal performance in the self-organizing model? It is likely to do with a lack of specificity in the COMB cell layer. The projections from a given COMB cell onto the HD layer will only be correct if that COMB cell is active in response to a single HD. COMB cells that are active for more than one specific HD can, for at least one of these HDs, introduce error into path integration. This leads to the prediction of a specific source of error within the HD system. However, it should be noted that, theoretically, COMB cell activity at multiple HDs would not necessarily be predicted to introduce error into the HD system: it would only be when COMB cell projections are close and strong enough to the current HD layer activity packet that they would affect packet motion.

The accuracy of path integration is an interesting question, not least because it is not clear exactly how accurate path integration is in the brain. One caveat in working out path integration error in the brain is that it tends to vary experimentally. Crucially, experimental measures of path integration error are linked to changes in allocentric preferred firing direction of a given HD cell. In the light, such changes would result from an altered relationship between HD activity and the external world resulting in remapped allothetic input (see, for example (Knight et al. 2014; Page et al. 2014)). It is therefore worth remembering that inaccuracy in HD packet motion can only be directly measured by shifts in preferred firing direction in the absence of allothetic information to which the HD system can orient.

In addition to potential sources of inaccuracy within the HD system (outlined here and in (Page et al. 2015)), there is a possibility of error being introduced prior to the HD system itself. It has been shown that rats tend to structure their exploratory movements in a highly specific manner. They start at a home base and perform a number of movements before returning home in a more direct fashion. During their outward journey, they typically segment their movements into paths of low curvature at higher speed, and paths of high curvature at a lower speed (i.e., the faster the rat travels, the straighter the path they take). This inverse correlation of speed and curvature likely enhances the gain of angular accelerations and makes them more detectable to the vestibular system (Wallace et al. 2006). For example, slower changes in heading whilst moving across a long curved path might be harder to detect than splitting the same path into a series of more abrupt changes with straighter paths. It has been shown in humans that the disruption of this inverse relationship is associated with errors in path integration (Wallace et al. 2006).

It could well be that movement segmentation falls apart due to errors in path integration, rather than causing them. However, it has been shown that when rats are forced to segment their path in a manner that disrupts inverse correlations of speed and path curvature, they are less accurate at returning to home than when they are free to structure their own movements (Wallace et al. 2010). It is likely that intrinsic HD system error can accumulate and cause a breakdown in the movement segmentation described above. This loss of movement segmentation compounds the error by reducing the accuracy of vestibular input. This would generate “upstream” errors in AHV input to the HD system.

Path integration accuracy is not yet fully understood from a computational perspective. It is non-trivial that errors in computational path integration systems tend to be undershooting (Walters et al. 2013; Page et al. 2015), yet experimental studies indicated preferred firing direction shifts in both directions (Goodridge et al. 1998); indicating that path integration errors are randomly distributed. Future studies of path integration ought to replicate this noisy element and move away from systematical undershooting of the desired speed. The understanding of path integration error resulting from a noisy biological process ought to be furthered, and for this reason, future work should focus on replicating experiments as closely as possible in silico. A target model framework to investigate path integration inaccuracy should incorporate a fully spiking network which simulates ecologically accurate movements at multiple velocities.

The simulations presented in this paper have been performed using a simplified rate-coded neural network model, in which the activity of individual neurons at any moment in time is represented as a mean firing rate. However, real neurons in the brain communicate with each other using discrete electrical pulses called action potentials or ‘spikes’. In future work, we plan to run simulations of the head direction system using a spiking neural network model that can more accurately mimic the temporal dynamics of neurons and synapses (Eguchi et al. 2018; Isbister et al. 2018). In particular, this will enable us to implement realistic neuronal and synaptic time constants governing, for example, the exponential decays of both cell voltages and synaptic conductances. These two-time constants will actually have different values in the brain. However, in the more simplified rate-coded simulations presented above, these two biological time constants are effectively combined into a single neuronal time constant thus reducing the temporal accuracy of the dynamics. In earlier work with spiking neural networks, we have implemented neuronal time constants of 20 ms and synaptic time constants of 1 ms. While in the rate coded simulations presented here, we explored network performance across a range of values for the neuronal time constants from 0.1 to 100 ms. So with the rate-coded simulations, the conclusions are more qualitative than quantitative. Nevertheless, despite the reduced dynamical accuracy of rate-coded models, the simulation results presented in this paper have demonstrated how removing the recurrent connections between head direction cells enables the model to produce arbitrarily accurate path integration, which now depends on the duration of the time constants. Shorter neuronal time constants produce faster neuronal rise (response) times, which in turn lead to more accurate path integration. In particular, the path integration error can be made arbitrarily small by reducing the neuronal time constants to very small values. We tested this observation by running simulations with a neuronal time constant reduced to an extremely low value of 0.1 ms. This finding represents an advance on previous continuous attractor models of the head direction system that relied on the recurrent connections between head direction cells to stabilize the activity packet in the absence of visual input (Walters and Stringer 2010; Walters et al. 2013). These earlier models suffered from much larger path integration inaccuracies caused by the recurrent connectivity, which could not be eliminated by adjusting model parameters such as the neuronal time constants. In particular, and most importantly, the simulation results presented in this paper are consistent with and provide an explanation for why the head direction system in the rat brain does not appear to incorporate recurrent connections between head direction cells. The reason that recurrent connections between head direction cells are absent in the rat brain is to help optimize the accuracy of path integration using the mechanisms described above in the Conclusions section.

Based on the known physiology of the HD cell system, the pre-wired and self-organizing models in this paper should be considered in terms of DTN and LMN, and thus we must attempt to match the connectivity of these regions. Roughly speaking, the HD region could be mapped onto the LMN, whilst the COMB cell region could be mapped onto the DTN. Unlike the connectivity in each of our models, which is uniformly excitatory, DTN to LMN connections appears to be inhibitory. In order to have strict biological accuracy, the models should, therefore, be extended to cover inhibitory connectivity from COMB cells back onto the HD layer. One potential issue with this change, particularly for self-organization, is how to update weights of inhibitory connections. One solution to this issue is the concept of anti-Hebbian learning, whereby weights are increased based on anticorrelation between pre- and post-synaptic firing rather than correlation (Koch et al. 2013; Vogels et al. 2013). However, as stated above, in the model presented in this paper the connections from the COMB layer to the HD layer were excitatory. In future work, we plan to investigate how the model may operate when these connections are inhibitory. That is, can inhibitory connections from the COMB layer to the HD layer self-organize using anti-Hebbian learning to allow the network to learn to perform accurate path integration?

The assertion in this model that it is the recurrent connectivity between DTN and LMN which constitutes the core HD system attractor could be tested by specifically disrupting reciprocal connectivity between these layers without damaging their internal connectivity. This would be analogous to severing the HD to COMB and COMB to HD connections within this model and represents a good test of this model’s plausibility. The other assertion in this model that HD signal generation is critically dependent on AHV input, has already been approached. (Valerio and Taube 2016) have shown that a stable HD signal is no longer present in the ADN of mice without AHV input, eliminated at the level of the semicircular canals. Instead, HD cells fire in bursts (i.e., with a pattern of spiking similar to HD cells, but not tied to a specific allocentric direction). This seems at odds with the current model, which suggests a loss of AHV input will silence the HD system. However, it is perfectly plausible that other, less primary, sources of idiothetic (e.g., optic flow or motor efference copy – see (Zugaro et al. 2001; Stackman et al. 2003)) and allothetic information (e.g., olfaction or audition – see (Goodridge et al. 1998)) contribute to HD firing, yet result in spatial instability of HD cell firing without AHV input. For example, if visual landmarks are not stably learnt with respect to idiothetic information then they will not have a strong or consistent influence on HD layer activity (Page et al. 2014). Therefore, it is likely that bursty firing arises from a combination of allocentric information sources which cannot be stably learnt, in conflict with idiothetic cues which themselves are not enough to stably influence HD packet motion.

The models presented here have attempted to simplify operation enough to gain insight without being lost in complexity. In particular, the firing of ROT cells is greatly simplified when compared to AHV cells in the brain. AHV cells have a roughly linear firing rate response across a range of head velocities and vary with their response to head rotation of different directions. Here, ROT cells fire in a binarized manner, either firing at maximum for the preferred velocity of head rotation or at zero otherwise. A natural extension of this work would, therefore, be to have ROT cells more accurately reflecting AHV cell firing. Given a more biologically accurate ROT cell layer with a less caricatured response to head rotational velocity, it would be interesting to see whether this network could still function. It would also be of interest to see what influence a population encoding in the ROT cell layer would have on the representations formed in the COMB cell layer.

It is notable that the COMB cells in this paper are also simplified relative to biology, in that they respond to a single HD and a single speed of rotation. In the brain, cells responding to combinations of HD and angular head velocity show elevated firing across a range of head directions, and this is modulated in a complex manner by head rotation (Bassett and Taube 2001). It would be of interest then, to see if realistic COMB cell representations follow from realistic ROT cell representations. The current models deal only with two rotation speeds, with one being zero, and with the non-zero speed being in a single direction. It is possible that a richer ROT/COMB cell representation covering different speeds would allow the system to function across a range of angular head velocities, and it would be interesting to see how the system copes with this increased challenge, as well as what sort of weight profiles would be required.

In a previous theoretical study (Page et al. 2015), we carried out a detailed simulation analysis of how symmetric recurrent connections, as set-up by a biologically plausible Hebbian associative learning rule, within an HD layer will lead to path integration error. These symmetric recurrent connections always slow the movement of the activity packet within the HD layer when the simulated agent rotates in the dark, thereby producing an underestimate of the angular change in head direction. More generally, the underlying problem is that in order to avoid path integration error during head rotation, the recurrent connections would need to be suitably shaped, e.g., with an appropriate offset around the HD ring. It is not possible to achieve this for a broad range of different rotation speeds simultaneously using the same set of recurrent connections, unless perhaps by applying some biologically implausible modulatory mechanism to these connections based on rotation speed. This theoretical observation would appear to offer an explanation for why there do not appear to be excitatory recurrent connections within individual layers of the head direction system in the rat brain. However, the model simulations in that previous study were performed using a simplified one-layer neural network. In the current paper, we have extended this work to a two-layer network and shown how the connections may self-organize using biologically plausible ‘local’ learning rules.

To conclude, here we present a two-layer model capable of sustaining a stationary activity packet, and shifting an activity packet using an internal head rotation signal with a high degree of accuracy. This path integration accuracy was achieved by eliminating the recurrent excitatory connections within the HD layer, which are present in many HD system models. Continuing sources of path integration inaccuracy are hypothesized: 1) neuronal rise time; 2) incorrect COMB layer influence on HD firing, due either to nonspecific HD influence on COMB, COMB influence on HD, or both; 3) ”upstream” errors in AHV input. This work represents a useful contribution to our theoretical understanding of the rat HD cell system, and in particular, provides a platform from which to investigate the self-organization of path integration and potential sources of error in more detail.

Acknowledgments

This work was funded by the Oxford Foundation for Theoretical Neuroscience and Artificial Intelligence.

Disclosure statement

No potential conflict of interest was reported by the authors.