Learning accurate path integration in ring attractor models of the head direction system

The gross model architecture closely resembles the one found in the fly CX (Figure 1A). It comprises HD cells organized in a ring, and HR cells organized in two wings. One wing is responsible for leftward and the other for rightward movement of the internal heading representation. HD cells receive visual input from the so-called ‘ring’ neurons; this input takes the form of a disinhibitory bump centered at the current HD (Figure 1B, Omoto et al., 2017; Fisher et al., 2019). The location of this visual bump in the network is controlled by the current head direction. We simulate head movements by sampling head-turning velocities from an Ornstein-Uhlenbeck process (Materials and methods), and we provide the corresponding velocity input to the HR cells (Figure 1C). HR cells provide direct input to HD cells, and HR cells also receive input from HD cells (Figure 1A). Both HR and HD cells receive global inhibition, which is in line with a putative ‘local’ model of HD network organization (Kim et al., 2017). The connections from HR to HD cells ($W^{\mathrm{HR}}$) and the recurrent connections among HD cells ($W^{\text{rec}}$) are assumed to be plastic. The goal of learning is to tune these plastic connections so that the network can achieve PI in the absence of visual input.

Figure 1 with 2 supplements see all

Download asset Open asset

The unit that controls plasticity in our network is an ‘associative neuron’. It is inspired by pyramidal neurons of the mammalian cortex whose dendrites act, via backpropagating action potentials, as coincidence detectors for signals arriving from different layers of the cortex and targeting different compartments of the neuron (Larkum et al., 1999). Paired with synaptic plasticity, coincidence detection can lead to long-lasting associations between these signals (Larkum, 2013). To map the morphology of a cortical pyramidal cell to the one of a HD cell in the fly, we first point out that all relevant inputs arrive at the dendrites of HD cells within the ellipsoid body (EB) of the fly (Xu, 2020) moreover, the soma itself is externalized in the fly brain, and it is unlikely to contribute considerably to computations (Gouwens and Wilson, 2009; Tuthill, 2009). We thus link the dendrites of the pyramidal associative neuron to the axon-distal dendritic compartment of the associative HD neuron in the fly, and we link the soma of the pyramidal associative neuron to the axon-proximal dendritic compartment of the associative HD neuron in the fly. Furthermore, we assume that the axon-proximal compartment is electrotonically closer to the axon initial segment, and therefore, similarly to the somatic compartment in pyramidal neurons, inputs there can more readily initiate action potentials. Note that our model does not require active backpropagation of action potentials — passive spread of voltage to the axon-distal compartment would be sufficient (for details, see Materials and methods and Discussion). We also assume that associative HD cells receive visual input ($I^{\text{vis}}$) in the axon-proximal compartment, and both recurrent input ($W^{\text{rec}}$) and HR input ($W^{\text{HR}}$) in the axon-distal compartment; accordingly, we model HD neurons as two-compartment units (Figure 1D). The associative neuron can learn the synaptic weights of the incoming connections in the axon-distal compartment, therefore, as mentioned, we let $W^{\text{rec}}$ and $W^{\text{HR}}$ be plastic.

We find that the assumption of spatial segregation of postsynapses of HD cells is consistent with our analysis of EM data from the fly (Xu, 2020). For an example HD (E-PG) neuron, Figure 1E depicts that head rotation and recurrent inputs (mediated by P-EN1 and P-EN2 cells, respectively [Turner-Evans et al., 2020]) contact the E-PG cell in locations within the EB that are distinct compared to those of visually responsive neurons R2 and R4d (Omoto et al., 2017; Fisher et al., 2019), as hypothesized. The same pattern was observed for a total of 16 E-PG neurons (one for each ‘wedge’ of the EB) that we analyzed (Figure 1—figure supplement 1A). To further support the assumption that visual inputs are separated from recurrent and HR-to-HD inputs, we perform binary classification between the two classes, using SVMs (for details, see Materials and methods). Figure 1—figure supplement 1B shows that predicting class identity from spatial location alone in held-out test data is excellent (test accuracy >0.95 across neurons and model runs).

The connections from HD to HR cells ($W^{\text{HD}}$) are assumed to be fixed, and HR cells are modeled as single-compartment units. Projections are organized such that each wing neuron receives input from only one specific HD neuron for every HD (Figure 1A). This simple initial wiring makes HR cells conjunctively tuned to HR and HD, and we assume that it has already been formed, for example, during circuit assembly. We note that the conditions for 1-to-1 wiring and constant amplitude of the HD-to-HR connections can be relaxed, because the learning rule can balance asymmetries in the initial architecture (see Appendix 3). In addition, the connections carrying the visual and angular velocity inputs are also assumed to be fixed. Although plasticity in the visual inputs has been shown to exist (Fisher et al., 2019; Kim et al., 2019), here we focus on how the path-integrating circuit itself originally self-organizes. Therefore, to simplify the setting and without loss of generality, we assume a fixed anchoring to environmental cues as the animal moves in the same environment (for details, see Discussion).

In our model, the visual input acts as a supervisory signal during learning (as in D’Albis and Kempter, 2020), which is used to change weights of synapses onto the axon-distal compartment of HD cells. We utilize the learning rule proposed by Urbanczik and Senn, 2014 (for details, see Materials and methods), which tunes the incoming synaptic connections in the axon-distal compartment in order to minimize the discrepancy between the firing rate of the neuron $f(V^a)$ (where $V^a$ is the axon-proximal voltage, primarily controlled by the visual input) and the prediction of the firing rate by the axon-distal compartment from axon-distal inputs alone, $f(p{V}^d)$ (where $p$ is a constant and $V^d$ is the axon-distal voltage, which depends on head rotation velocity). From now on, we refer to this discrepancy as ‘learning error’, or simply ‘error’ (Equation 18; in units of firing rate). The synaptic weight change $\mathrm{\Delta }{W}_{\mathit{\text{pre,post}}}$ from a presynaptic (HD or HR) neuron to a postsynaptic HD neuron is then given by:

where $\eta$ is the constant learning rate and ${\text{P}}_{\text{pre}}$ is the postsynaptic potential from the presynaptic neuron. When implementing this learning rule, we low-pass filter the prospective weight change $\mathrm{\Delta }{W}_{\mathit{\text{pre,post}}}$ to ensure smoothness of learning.

Importantly, this learning rule is biologically plausible because the firing rate of an associative neuron $f(V^a)$ is locally available at every synapse in the axon-distal compartment due to the (passive or active) backpropagation of axonal activity to the axon-distal dendrites. The other two signals that enter the learning rule are the voltage of the axon-distal compartment $V^d$ and the postsynaptic potential P, which are also available locally at the synapse; for details, see Materials and methods. Furthermore, recent behavioral experiments show that conditioning in Drosophila (Zhao et al., 2021) is not well explained by classical correlation-based plasticity, but it can be well accounted for by predictive synaptic plasticity. The latter is in line with the learning rule utilized here.

Figure 2A shows an example of the performance of a trained network, for the light condition (i.e. when visual input is available; yellow overbars) and for PI in darkness (purple overbars); the performance is quantified by the PI error (in units of degrees) over time. PI error refers to the accumulated difference between the internal representation of heading and the true heading, and it is different from the learning error introduced previously.

Figure 2

Download asset Open asset

A unique bump of activity is clearly present at all times in the HD network (Figure 2A, top), in both light and darkness conditions, and this bump moves smoothly across the network for a variable angular velocity (Figure 2A, bottom). The position of the bump is defined as the population vector average (PVA) of the neural activity in the HD network. The HD bump also leads to the emergence of bumps in the HR network, separately for L-HR and R-HR cells (Figure 2A, second and third panel from top). In light conditions (0–20 s in Figure 2A), the PVA closely tracks the head direction of the animal in HD, L-HR, and R-HR cells alike, which is expected because the visual input guides the network activity. Importantly, however, in darkness (20–50 s in Figure 2A), the self-motion input alone is enough to track the animal’s heading, leading to a small PI error between the internal representation of heading and the ground truth. This error is corrected after the visual input reappears (at 50 s in Figure 2A). Such PI errors in darkness are qualitatively consistent with data reported in the experimental literature (Seelig and Jayaraman, 2015). The correction of the PI error also reproduces in silico the experimental finding that the visual input (whenever available) exerts stronger control on the bump location than the self-motion input (Seelig and Jayaraman, 2015), which suggests that even the mature network does not rely on PI when visual cues are available.

To quantify the accuracy of PI in our model, we draw 1,000 trials, each 60 s long, for constant synaptic weights and in the absence of visual input. We also limit the angular velocities in these trials to retain only velocities that flies realistically display (see dashed green lines in Figure 2C and Methods). We then plot the distribution of PI errors every 10 s (Figure 2B). We find that average absolute PI errors (widths of distributions) increase with time in darkness, but most of the PI errors at 60 s are within 60 deg of the true heading. This vastly exceeds the PI performance of flies (Seelig and Jayaraman, 2015). In flies, the correlation between the PVA estimate and the true heading in darkness varied widely across animals in the range [0.3, 0.95] (Seelig and Jayaraman, 2015), whereas for the model it is close to 1. However, it should be noted that the model here corresponds to an ideal scenario that serves as a proof of principle. We will later incorporate irregularities owing to biological factors (asymmetry in the weights, biological noise) that bring the network’s performance closer to the fly’s behavior.

To further assess the network’s ability to integrate different angular velocities, we simulate the system both with and without visual input in 5 s intervals during which the angular velocity is constant. We then compute the average movement velocity of the bump across the network, that is the neural velocity, and compare it to the real velocity provided as input. Figure 2C shows that the network achieves a PI gain (defined as the ratio between neural and real velocity) close to 1 both with and without supervisory visual input, meaning that the neural velocity matches very well the angular velocity of the animal, for all angular velocities that are observed in experiments (${\displaystyle |v|<500}$ deg/s for walking and flying) (Geurten et al., 2014; Stowers et al., 2017). Although expected in light conditions, the fact that gain 1 is achieved in darkness shows that the network predicts the missing visual input from the velocity input, that is, the network path integrates accurately. Note that PI is impaired in our model for very small angular velocities (Figure 2C, flat purple line for ${\displaystyle |v|<30}$ deg/s), similarly to previous hand-tuned theoretical models (Turner-Evans et al., 2017). This is a direct consequence of the fact that maintaining a stable activity bump and moving it across the network at very small angular velocities are competing goals. Crucially, it has been reported that such an impairment of PI for small angular velocities exists in flies (Seelig and Jayaraman, 2015). Note that if we increase the number of HD neurons from 60 (∼50 were reported in the fly by Turner-Evans et al., 2020; Xu, 2020) to 120 or 240, this flat region is no longer observed (data not shown).

A continuous attractor network (CAN) should be able to maintain a localised bump of activity in virtually a continuum of locations around the ring of HD cells. To prove that the learned network approximates this property, we seek to reproduce in silico experimental findings in Kim et al., 2017. There it was shown that local optogenetic stimulation of HD cells in the ring can cause the activity bump to jump to a new position and persist in that location — supported by internal dynamics alone.

To reproduce the experiments by Kim et al., 2017, we simulate optogenetic stimulation of HD cells in our network as visual input of increased strength and extent (for details, see Materials and methods). We find that the strength and extent of the stimulation needs to be increased relative to that of the visual input; only in this case, a bump at some other location in the network can be suppressed, and a new bump emerges at the stimulated location. The stimuli are assumed to appear instantaneously at random locations, but we restrict our set of stimulation locations to the discrete angles represented by the finite number of HD neurons. Furthermore, the velocity input is set to zero for the entire simulation, signaling lack of head movement.

Figure 2D shows network activity in response to several stimuli, when the stimulation location changes abruptly every 5 s. During stimulation (2 s long, red overbars), the bump is larger than normal due to the use of a stronger than usual visual-like input to mimic optogenetic stimulation. The way in which the network responds to a stimulation depends on how far away from the ‘current’ location it is stimulated: for shorter distances, the bump activity shifts to the new location, as evidenced by the transient dynamics at the edges of the bump resembling a decay from an initial to a new location (see Figure 2D at {5,15,20} s). However, for larger phase shifts $\mathrm{\Delta }\theta$ the bump first emerges in the new location and subsequently disappears at the initial location, a mechanism akin to a ‘jump’ (Figure 2D, all other transitions). Similar effects have been observed in the experimental literature (Seelig and Jayaraman, 2015; Kim et al., 2017). The way the network responds to stimulation indicates that it operates in a CAN manner, and not as a winner-takes-all network where changes in bump location would always be instantaneous (Carpenter and Grossberg, 1987; Itti et al., 1998; Wang, 2002). That is to say, the network operates as expected from a quasi-continuous attractor. Furthermore, we find that the transition strategy in our model changes from predominantly smooth transitions to jumps at $\mathrm{\Delta }\theta \approx 90\text{deg}$, which matches experiments well (Kim et al., 2017).

Following a 2 s stimulation, the network activity has converged to the new cued location. After the stimulation has been turned off, the bump remains at the new location (within the angular resolution $\mathrm{\Delta }\varphi$ of the network), supported by internal network dynamics alone (Figure 2D). We confirmed in additional simulations that the bump does not drift away from the stimulated location for extended periods of time (3 min duration tested, only 3 s shown), and for all discrete locations in the HD network (only six locations shown). Therefore, we conclude that the HD network is a quasi-continuous attractor that can reliably sustain a heading representation over time in all HD locations. Note that for the network size used ($N^{\text{HD}}=60$) we still obtain discrete attractors with separated basins of attraction; however it is expected that with increasing $N^{\text{HD}}$ adjacent attractors will merge when the intrinsic noise overcomes the barrier separating them. Indeed, we find that for $N^{\text{HD}}=N^{\text{HR}}=120$ it is easier to diffuse to adjacent attractors in the presence of synaptic input noise; for the impact of noise, see Appendix 1—figure 1C. In reality, the bump may drift away due to asymmetries in the connectivity of the biological circuit as well as intrinsic noise (Burak and Fiete, 2012) see also Appendix 1. In flies, for instance, the bump can stay put only for several seconds (Kim et al., 2017).

To gain more insight into how the network achieves PI and attains CAN properties, we show how the synaptic weights of the network are tuned during a developmental period (Figure 3). Figure 3A and B shows the learned recurrent synaptic weights among the HD cells, $W^{\text{rec}}$, and the learned synaptic weights from HR to HD cells, $W^{\text{HR}}$, respectively. Circular symmetry is apparent in both matrices, a crucial property for a symmetric ring attractor. Therefore, we also plot the profiles of the learned weights as a function of receptive field difference in Figure 3C. Note that pixelized appearance in these plots is due to the fact that two adjacent HD neurons are tuned for the same HD, and develop identical synaptic strengths.

Figure 3 with 3 supplements see all

Download asset Open asset

First, we discuss the properties of the learned weights. Local excitatory connections have developed along the main diagonal of $W^{\text{rec}}$, similar to what is observed in the CX (Turner-Evans et al., 2020). This local excitation can be readily seen in the weight profile of $W^{\text{rec}}$ in Figure 3C, and it is the substrate that allows the network to support stable activity bumps in virtually any location. In addition, we observe inhibition surrounding the local excitatory profile in both directions. This inhibition emerges despite the fact that we provide global inhibition to all HD cells ($I_{inh}^{\text{HD}}$ parameter, Materials and methods), in line with suggestions from previous work (Kim et al., 2017). Surrounding inhibition was a feature we observed consistently in learned networks of different sizes and for different global inhibition levels. Finally, the angular offset of the two negative sidelobes in the connectivity depends on the size and shape of the entrained HD bump (for details, see Appendix 5).

Furthermore, we find a consistent pattern of both L-HR and R-HR populations to excite the direction for which they are selective (Figure 3C), which is also similar to what is observed in the CX (Turner-Evans et al., 2020). Excitation in one direction is accompanied by inhibition in the reverse direction in the learned network. As a result of the symmetry in our learning paradigm, the connectivity profiles of L-HR and R-HR cells are mirrored versions of each other, which is also clearly visible in Figure 3C. The inhibition of the reverse direction has a width comparable to the bump size and acts as a ‘break’ to prevent the bump from moving in this direction. The excitation in the selective direction, on the other hand, has a wider profile, which allows the network to path integrate for a wide range of angular velocities, that is for high angular velocities neurons further downstream can be ‘primed’ and activated in rapid succession. Indeed, when we remove the wide projections from the excitatory connectivity, PI performance is impaired for the higher angular velocities exclusively (Figure 3—figure supplement 1). The even weight profile in $W^{\text{rec}}$ and the mirror symmetry for L-HR vs. R-HR profiles in $W^{\text{HR}}$, together with the circular symmetry of the weights throughout the ring, guarantee that there is no side bias (i.e. tendency of the bump to favor one direction of movement versus the other) during PI. Indeed, the PI error distribution in Figure 2B remains symmetric throughout the 60 s simulations.

Next, we focus our attention on the dynamics of learning. For training times larger than a few hours, the absolute learning error drops and settles to a low value, indicating that learning has converged after ∼20 hr (or 4000 cycles, each cycle lasting $1/\eta$) of training time (Figure 3D). The non-zero value of the final error is only due to errors occurring at the edges of the bump (Figure 3—figure supplement 2A, top panel). An intuitive explanation of why these errors persist is that the velocity pathway is learning to predict the visual input; as a result, when the visual input is present, the velocity pathway creates errors that are consistent with PI velocity biases in darkness.

Figure 3E and F shows the weight development history for the entire simulation. The first structure that emerges during learning is the local excitatory recurrent connections in $W^{\text{rec}}$. For these early stages of learning, the initial connectivity is controlled by the autocorrelation of the visual input, which gets imprinted in the recurrent connections by means of Hebbian co-activation of adjacent HD neurons. As a result, the width of the local excitatory profile mirrors the width of the visual input. Once a clear bump is established in the HD ring, the HR connections are learned to support bump movement, and negative sidelobes in $W^{\text{rec}}$ emerge. To understand the shape of the learned connectivity profiles and the dynamics of their development, we study a reduced version of the full model, which follows learning in bump-centric coordinates (see Appendix 5). The reduced model produces a connectivity strikingly similar to the full model, and highlights the important role of non-linearities in the system.

So far, we have shown results in which our model far outperforms flies in terms of PI accuracy. To bridge this gap, we add noise to the weight connectivity in Figure 3A and B and obtain the connectivity matrices in Figure 3—figure supplement 3A,B, respectively. This perturbation of the weights could account for irregularities in the fly HD system owning to biological factors such as uneven synaptic densities. The resulting neural velocity gain curve in Figure 3—figure supplement 3E is impaired mainly for small angular velocities (Figure 2C). Interestingly, it now bears greater similarity to the one observed in flies, because the previously flat area for small angular velocities is wider (flat for ${\displaystyle |v|<60}$ deg/s, cf. extended data fig. 7G,J in Seelig and Jayaraman, 2015). This happens because the noisy connectivity is less effective in initiating bump movement. Finally, the PI errors in the network with noisy connectivity grow much faster and display a strong side bias (Figure 3—figure supplement 3D, Figure 2B). The latter can be attributed to the fact that the noise in the connectivity generates local minima that are easier to transverse from one direction vs. the other. Side bias can also emerge if the learning rate $\eta$ in Equation 16 is increased, effectively forcing learning to converge faster to a local minimum, which results in slight deviations from circularly symmetric connectivity (data not shown). It is therefore expected that different animals will display different degrees and directions of side bias during PI, owning either to fast learning or asymmetries in the underlying neurobiology. Since the exact behavior of the network with noise in the connectivity depends on the specific realization, we also generate multiple such networks and estimate the diffusion coefficient during path integration, which quantifies how fast the width of the PI error distribution in Figure 3—figure supplement 3D increases. We find the grand average to be $82.3\pm 15.7{\text{deg}}^2/\text{s}$, which is considerably larger (Student’s t-test, 95% conf. intervals for a total of 12 networks) than the diffusion coefficient for networks without a perturbation in the weights ($24.5{\text{deg}}^2/\text{s}$ in Figure 2B). Finally, in Appendix 1 we also incorporate random Gaussian noise to all inputs, which can account for noisy percepts or stochasticity of spiking, and show that learning is not disrupted even for high noise levels.

Having shown how PI and CAN properties are learned in our model, we now turn our attention to the flexibility that our learning setup affords. Motivated by augmented-reality experiments in rodents where the relative gain of visual and self-motion inputs is manipulated (Jayakumar et al., 2019), we test whether our network can rewire to learn an arbitrary gain between the two. In other words, we attempt to learn an arbitrary gain $g$ between the idiothetic angular velocity $v$ sensed by the HR cells and the neural velocity $g\cdot v$ dictated by the allothetic visual input. This simulates the conditions in an augmented reality environment, where the speed at which the world around the animal rotates is determined by the experimenter, but the proprioceptive sense of head angular velocity remains the same.

Starting with the learned network shown in Figure 3, which displayed gain $g=1$, we suddenly switch to a different gain, that is we learn weights for $g\in \{0.25,0.5,1.5,2\}$. In all cases, we observe that the network readily rewires to achieve the new gain. The mean learning error after the gain switch is initially high, but reaches a lower, constant level after at most 3 hr of training (Figure 4A). We note that convergence is much faster compared to the time it takes for the gain-1 network to emerge from scratch (compare to Figure 3D), especially for the smaller gain changes. Importantly, Figure 4B shows that PI performance in the resulting networks is excellent for the new gains, with some degradation only for very low and very high angular velocities. There are two reasons why high angular velocities are not learned that well: limited training of these velocities, and saturation of HR cell activity. Both reasons are by design and do not reflect a fundamental limit of the network.

Figure 4 with 1 supplement see all

Download asset Open asset

In Appendix 2, we show that without the aforementioned limitations the network learns to path-integrate up to an angular velocity limit set by synaptic delays and that the bump width sets a trade-off between location and velocity-integration accuracy in the HD system.

Figure 4C and D compare the weight profiles of the circularly symmetric matrices $W^{\text{rec}}$ and $W^{\text{HR}}$ resulting from the initial gain $g=1$, with the weight profiles resulting from adaptation to the most extreme gains shown in Figure 4, that is $g\in \{0.25,2\}$. An increase in gain slightly suppresses the recurrent connections and slightly amplifies the HR-to-HD connections, while a decrease in gain substantially amplifies the recurrent connections and slightly suppresses the HR-to-HD connections. The latter explains why the flat region for small angular velocities in Figure 4B has been extended for $g\in \{0.25,0.5\}$: it is now harder for small angular velocities to overcome the attractor formed by stronger recurrent weights and move the bump.

Finally, we address the limits of the ability of the network to rewire to new gains (Figure 4—figure supplement 1). We find that after rewiring the performance is excellent for gains between 0.25 and 4.5. The network can even reverse its gain to $g=-1$, that is, when allothetic and idiothetic inputs are signaling movement in opposite directions. However, for larger gain changes, learning takes longer.

Inspired by the recent discovery of a ring attractor network for HD in Drosophila (Seelig and Jayaraman, 2015), we show how a biologically plausible learning rule leads to the emergence of a circuit that achieves gain-1 PI in darkness. The learned network features striking similarities in terms of connectivity to the one experimentally observed in the fly (Turner-Evans et al., 2020), and reproduces experiments on CAN dynamics (Kim et al., 2017) and gain changes between external and self-motion cues in rodents (Jayakumar et al., 2019). Furthermore, an impairment of PI for small angular velocities is observed in the mature network, which is a feature that has been reported in experiments (Seelig and Jayaraman, 2015). Finally, the proposed learning rule can serve to compensate deviations from circular symmetry in the synaptic weight profiles; such deviations are expected in biological systems and — if not compensated — could lead to large PI errors.

The mature circuit displays two properties characteristic of CANs: (1) it can support and actively maintain a local bump of activity at a virtual continuum of locations, and (2) it can move the bump across the network by integrating self-motion cues. Note that we did not explicitly train the network to achieve these CAN properties, but they rather emerged in a self-organized manner.

To achieve gain-1 PI performance, our network must attribute learning errors to the appropriate weights. The learning rule we adopt in Equation 1 is a ‘delta-like’ rule, with a learning error that gates learning in the network, and a Hebbian component that comes in the form of the postsynaptic potential and assigns credit to synapses that are active when errors are large. The learning rule leads to the emergence of both symmetric local connectivity between HD cells (which is required for bump maintenance and stability), and asymmetric connectivity from HR to HD cells (which is required for bump movement in darkness). The first happens because adjacent neurons are co-active due to correlated visual input; the second because only one HR population is predominantly active during rotation: the population that corresponds to the current rotation direction. Crucial to the understanding of the learning dynamics of the model was the development of a reduced model, which follows learning in bump-centric coordinates and is analytically tractable (see Appendix 5). The reduced model can be extended to higher dimensional manifolds (Gardner et al., 2022), and therefore it offers a general framework to study how activity-dependent synaptic plasticity shapes CANs.

Our work comes at a time at which the fly HD system receives a lot of attention (Seelig and Jayaraman, 2015; Turner-Evans et al., 2020; Kim et al., 2017; Kim et al., 2019; Fisher et al., 2019), and suggests a mechanism of how this circuit could self-organize during development. Synaptic plasticity has been shown to be important in this circuit for anchoring the visual input to the HD neurons when the animal is exposed to a new environment (Kim et al., 2019; Fisher et al., 2019). This has also been demonstrated in models of the mammalian HD system (Skaggs et al., 1995; Zhang, 1996; Song and Wang, 2005). Here, we assume that an initial anchoring of the topographic visual input to the HD neurons with arbitrary offset with respect to external landmarks already exists prior to the development of the PI circuit; such an anchoring could even be prewired. In our model, it is sufficient that the visual input tuning is local and topographically arranged. Once the PI circuit has developed, visual connections could be anchored to different environments, as shown by Kim et al., 2019 and Fisher et al., 2019. Alternatively, the HD system itself could come prewired with an initial gross connectivity, sufficient to anchor the visual input; in this case, our learning rule would enable fine tuning of this connectivity for gain-1 PI. In either case, for the sake of simplicity and without loss of generality, we study the development of the path-integrating circuit while the animal moves in the same environment, and keep the visual input tuning fixed. Therefore, the present work addresses the important question of how the PI circuit itself could be formed, and it is complementary to the problem of how allothetic inputs to the PI circuit are wired (Fisher et al., 2019; Kim et al., 2019). The interplay of the two forms of plasticity during development would be of particular future interest.

A requirement for the learning rule we use is that information about the firing rate of HD neurons is available at the axon-distal compartment. There is no evidence for active backpropagation of APs in E-PG neurons in the fly, but passive backpropagation would suffice in this setting. In fact, passive spread of activity has been shown to attenuate weakly in central fly neurons (Gouwens and Wilson, 2009). In HD neurons, the axon-proximal and axon-distal compartments belong to the same dendritic tuft (Figure 1E), and since we assume that the axon initial segment is close to the axon-proximal compartment, the generated AP would need to propagate only a short distance compared to the effective electrotonic length. This means that APs would not be attenuated much on their way from the axon initial segment to the axon-distal compartment, and thus would maintain some of their high-frequency component, which could be used at synapses to differentiate them from slower postsynaptic potentials.

In Figure 4, we show that our network can adapt to altered gains much faster than the time required to learn the network from scratch. Our simulations are akin to experiments where rodents are placed in a VR environment and the relative gain between visual and proprioceptive signals is altered by the experimenter (Jayakumar et al., 2019). In this scenario, Jayakumar et al., 2019 found that the PI gain of place cells can be recalibrated rapidly. In contrast, Seelig and Jayaraman, 2015 found that PI gain in darkness is not significantly affected when flies are exposed to different gains in light conditions. We note, however, that Seelig and Jayaraman, 2015 tested mature animals (8–11 days old), whereas plasticity in the main HD network is presumably stronger in younger animals. Also note that the manipulation we use to address adaptation of PI to different gains differs from the one in Kim et al., 2019 who used optogenetic stimulation of the HD network combined with rotation of the visual scene to trigger a remapping of the visual input to the HD cells in a Hebbian manner. The findings in Jayakumar et al., 2019 can only be reconciled by plasticity in the PI circuit, and not in the sensory inputs to the circuit.

In order to address the core mechanisms that underlie the emergence of a path integrating network, we use a model that is a simplified version of the biological circuit. For example, we did not model inhibitory neurons explicitly and omitted some of the recurrent connectivity in the circuit, whose functional role is uncertain (Turner-Evans et al., 2020). We also choose to separate PI from other complex processes that occur in the CX (Raccuglia et al., 2019). Finally, we do not force the network to obey Dale’s law and do not model spiking explicitly.

Nevertheless, after learning, we obtain a network connectivity that is strikingly similar to the one of the fly HD system. Indeed, the mature model exhibits local excitatory connectivity in the HD neurons (Figure 3A and C), which in the fly is mediated by the excitatory loop from E-PG to P-EG to P-EN2 and back to E-PG (Turner-Evans et al., 2020), a feature that hand-tuned models of the fly HD system did not include (Turner-Evans et al., 2017). Furthermore, the HR neurons have excitatory projections towards the directions they are selective for (Figure 3B and C), similar to P-EN1 neurons in the fly. Interestingly, these key features that we uncover from learning have been utilized in other hand-tuned models of the system (Turner-Evans et al., 2017; Kim et al., 2017; Kim et al., 2019). Future work could endeavor to come closer to the architecture of the fly HD system and benefit from the incorporation of more neuron types and the richness of recurrent connectivity that has been discovered in the fly (Turner-Evans et al., 2020).

Compared to the fly, our network achieved better PI performance. As a simple way to match the performances, we added noise to the learned connectivity in the model; however, this is not an explanation why the fly performs worse. Indeed, there could be multiple reasons why PI performance is worse in the biological circuit. For instance, a confounder that would affect performance but not necessarily learning could be the presence of inputs that are unrelated to path integration, for example, inputs related to circadian cycles and sleep (Raccuglia et al., 2019). In the presence of such confounders, a precise tuning of the weights might be crucial in order to reach the performance of the fly. In other words, only if the model outperforms the biological circuit in a simplified setting, it has a chance to perform as well in a realistic setting, with all the additional complexities the latter comes with.

A common problem with CANs is that they require fine tuning: even a slight deviation from the optimal synaptic weight tuning leads to catastrophic drifting (Goldman et al., 2009). A way around this problem is to sacrifice the continuity of the attractor states in favor of a discrete number of stable states that are much more robust to noise or weight perturbations (Kilpatrick et al., 2013). In our network, the small number of HD neurons enables a coarse-grained representation of heading; the network is a CAN only in a quasi-continuous manner, and the number of discrete attractors corresponds to the number of HD neurons. This makes it harder to transition to adjacent attractors, since a ‘barrier’ has to be overcome in the quasi-continuous case (Kilpatrick et al., 2013). The somewhat counter-intuitive conclusion follows that a CAN with more neurons and, as a result, finer angular resolution, will not be as potent in maintaining activity, and diffusion to nearby attractors will be easier since the barrier will be lower. Indeed, we found that doubling the number of neurons produces a CAN that is less robust to noise. Overall, the quasi-continuous and coarse nature of the attractor shields the internal representation of heading against the ever-present biological noise, which would otherwise lead to diffusion of the bump with time. The fact that the network can still path-integrate accurately with this coarse-grained representation of heading is remarkable.

Seminal theoretical work on ring attractors has proven that in order to achieve gain-1 PI, the asymmetric component of the network connectivity (corresponding here to $W^{\text{HR}}$) needs to be proportional to the derivative of the symmetric component (corresponding to $W^{\text{rec}}$) (Zhang, 1996). However, this result rests on the assumption that asymmetric and symmetric weight profiles are mediated by the same neuronal population, as in the double-ring architecture proposed by Xie et al., 2002 and Hahnloser, 2003, but does not readily apply to the architecture of the fly HD system where HD and HR cells are separate. In our learned network, we find that the HR weight profile is not proportional to the derivative of the recurrent weight profile, therefore this requirement is not necessary for gain-1 PI in our setting. Note that our learning setup can also learn gain-1 PI for a double-ring architecture, which additionally obeys Dale’s law (Vafidis, 2019). Finally, we emphasize that circular symmetry is not a necessary condition for a ring attractor (Darshan and Rivkind, 2021). Rather, symmetry in our model results from the symmetry in the architecture, the symmetrically prewired weights, and the symmetric stimulus space. If any of those were to be relaxed, the resulting network would not be circular symmetric; then, the reduced model analysis that we perform in Appendix 5 would also not be feasible, because local asymmetries in the setup would result in non-local deviations from circular symmetry of the learned weights, which was our main assumption there. Nevertheless, we demonstrated that the full model can handle such asymmetries in the setup and learn accurate PI (see Appendix 3).

Our learning setup, inspired by Urbanczik and Senn, 2014, is similar to the one in Guerguiev et al., 2017 in the sense that both involve compartmentalized neurons that receive ‘target’ signals in a distinct compartment. It differs, however, in the algorithm and learning rule used. Guerguiev et al., 2017 use local gradient descent during a ‘target’ phase, which is separate from a forward propagation phase, akin to forward/backward propagation stages in conventional deep learning. In contrast, we use a modified Hebbian rule, and in our model ‘forward’ computation and learning happen at the same time; time multiplexing, whose origin in the brain is unclear, is not required. Our setting would be more akin to the one in Guerguiev et al., 2017 if an episode of PI in darkness would be required before an episode of learning in light conditions, which does not seem in line with the way animals naturally learn.

Previous theoretical work showed that head direction cells, head rotation cells, and grid cells emerge in neural networks trained for PI (Banino et al., 2018; Cueva and Wei, 2018). These networks were trained with backpropagation, therefore achieving gain-1 PI was not their primary focus; rather, this work elegantly demonstrated that the aforementioned cell types are efficient representations for spatial navigation that could be learned from experience.

We devote this section to discussing predictions of our model, and we suggest future experiments in flies and, potentially, other animal models. An obvious prediction of our model is that synaptic plasticity is critical for the development of the PI network for heading, and the lack of a supervisory allothetic sensory input (e.g. visual) during development should disrupt the formation of the PI system. Previous experimental work showed that head direction cells in rat pups displayed mature properties already in their first exploration of the environment outside their nest (Langston et al., 2010), which may seem to contradict our assumption that the PI circuit wires during development; however, directional selectivity of HD cells in the absence of allothetic inputs and PI performance were not tested in this study. In addition, it has been shown that visually impaired flies were not able to learn to accurately estimate the size of their body. This type of learning also requires visual inputs and, upon consolidation, remains stable (Krause et al., 2019).

We also predict that HD neurons have a compartmental structure where idiothetic inputs are separated from allothetic sensory inputs, which initiate action potentials more readily due to being electrotonically closer to the axon initial segment. While we already demonstrate the separation of allothetic and idiothetic inputs to E-PG neurons in the fly EB (Figure 1E, Figure 1—figure supplement 1), our prediction can only be tested with electrophysiological experiments. Another model prediction that can be tested only with electrophysiology is that APs backpropagate from the axon-proximal compartment (at least passively but with little attenuation) to the axon-distal compartment. Then spikes could be separated from postsynaptic potentials locally at the synapse by cellular mechanisms sensitive to the spectral density of the voltage.

Finally, similarly to place cell studies in rodents (Jayakumar et al., 2019), we predict that during development the PI system can adapt to experimenter-defined gain manipulations, and that it can do so faster than the time required for the system to develop from scratch. Therefore, a suggestion from this study would be to repeat in young flies the adaptation experiments by Seelig and Jayaraman, 2015.

The present study adds to the growing literature of potential computational abilities of compartmentalized neurons (Poirazi et al., 2003; Gidon et al., 2020; Payeur et al., 2021). The associative HD neuron used in this study is a coincidence detector, which serves to associate external and internal inputs arriving at different compartments of the cell. Coupled with memory-specific gating of internally generated inputs, coincidence detection has been suggested to be the fundamental mechanism that allows the mammalian cortex to form and update internal knowledge about external contingencies (Doron et al., 2020; Shin et al., 2021). This structured form of learning does not require engineered ‘hints’ during training, and it might be the reason why neural circuits evolved to be so efficient at reasoning about the world, with the mammalian cortex being the pinnacle of this achievement. Here, we demonstrate that learning at the cellular level can predict external inputs (visual information) by associating firing activity with internally generated signals (velocity inputs) during training. This effect is due to the anti-Hebbian component of the learning rule in Equation 12, where the product of postsynaptic axon-distal and presynaptic activity comes with a negative sign. Specifically, it has previously been demonstrated that anti-Hebbian synaptic plasticity can stabilize persistent activity (Xie and Seung, 2000) and perform predictive coding (Bell et al., 1997; Hahnloser, 2003). At the population level, this provides a powerful mechanism to internally produce activity patterns that are identical to the ones induced from an external stimulus. This mechanism can serve as a way to anticipate external events or, as in our case, as a way of ‘filling in’ missing information in the absence of external inputs.

Local, Hebb-like learning rules are considered a weak form of learning, due to their inability to utilize error information in a sophisticated manner. Despite that, we show that local associative learning can be particularly successful in learning appropriate fine-tuned synaptic connectivity, when operating within a cell structured for coincidence detection. Therefore, in learning and reasoning about the environment, our study highlights the importance of inductive biases with developmental origin (e.g. allothetic and idiothetic inputs arrive in different compartments of associative neurons) (Lake et al., 2017).

In conclusion, the present work addresses the age-old question of how to develop a CAN that performs accurate, gain-1 PI in the absence of external sensory cues. We show that this feat can be achieved in a network model of the HD system by means of a biologically plausible learning rule at the cellular level. Even though our network architecture is tailored to the one of the fly CX, the learning setup where idiothetic and allothetic cues are associated at the cellular level is general and can be applied to other PI circuits. Of particular interest is the rodent HD system: despite the lack of evidence for a topographically organized recurrent HD network in rodents, a one-dimensional HD manifold has been extracted in an unsupervised way (Chaudhuri et al., 2019). Therefore, our work lays the path to study the development of ring-like neural manifolds in mammals. Finally, it has recently been shown that grid cells in mammals form a continuous attractor manifold with toroidal topology (Gardner et al., 2022). It would be interesting to see if a similar mechanism underlies the emergence of PI in place and grid cells. Our model can be extended to higher dimensional CAN manifolds and provides a framework to interrogate this assumption.

We model a recurrent neural network comprising $N^{\text{HD}}=60$ head-direction (HD) and $N^{\text{HR}}=60$ head-rotation (HR) cells, which are close to the number of E-PG and P-EN1 cells in the fly central complex (CX), respectively (Turner-Evans et al., 2020; Xu, 2020). A scaled-down version of the network for $N^{\text{HR}}=N^{\text{HD}}=12$ is shown in Figure 1A. The average spiking activity of HD and HR cells is modelled by firing-rate neurons. HD cells are organized in a ring and receive visual input, which encodes the angular position of the animal’s head with respect to external landmarks. We use a discrete representation of angles and we model two HD cells for each head direction, as observed in the biological system (Turner-Evans et al., 2017). Therefore the network can represent head direction with an angular resolution $\mathrm{\Delta }\varphi =12\text{deg}$.

Motivated by the anatomy of the fly CX (Green et al., 2017; Turner-Evans et al., 2020), HR cells are divided in two populations (Figure 1A): a ‘leftward’ (L-HR) population (with increased velocity input when the head turns leftwards) and a ’rightward’ (R-HR) population (with increased velocity input when the head turns rightwards). After learning, these two HR populations are responsible to move the HD bump in the anticlockwise and clockwise directions, respectively.

The recurrent connections among HD cells and the connections from HR to HD cells are assumed to be plastic. On the contrary, connections from HD to HR cells are assumed fixed and determined as follows: for every head direction, one HD neuron projects to a cell in the L-HR population, and the other to a cell in the R-HR population. Because HD cells project to HR cells in a 1-to-1 manner, each HR neuron is simultaneously tuned to a particular head direction and a particular head rotation direction. The synaptic strength of the HD-to-HR projections is the same for all projections (these restrictions on the HD-to-HR connections are relaxed in Appendix 3). Finally, HR cells do not form recurrent connections.

We assume that each HD neuron is a rate-based associative neuron (Figure 1D), that is, a two-compartmental neuron comprising an axon-proximal and an axon-distal dendritic compartment (Urbanczik and Senn, 2014; Brea et al., 2016). The two compartments model the dendrites of that neuron that are closer to or further away from the axon initial segment. Note that here the axon-proximal compartment replaces the somatic compartment in the original model by Urbanczik and Senn, 2014. This is because the somata of fly neurons are typically electrotonically segregated from the rest of the cell and they are assumed to contribute little to computation (Gouwens and Wilson, 2009; Tuthill, 2009). We also note that to fully capture the input/output transformations that HD neurons in the fly perform, more compartments than two might be needed (Xu, 2020). Finally, only HD cells are associative neurons, whereas HR cells are simple rate-based point neurons.

HD cells receive an input current ${\mathit{I}}^d$ to the axon-distal dendrites, which obeys

where ${\mathit{I}}^d$ is a vector of length $N^{\text{HD}}$ with each entry corresponding to one HD cell. In Equation 2, ${\tau }_s$ is the synaptic time constant, $W^{\text{rec}}$ is a $N^{\text{HD}}\times N^{\text{HD}}$ matrix of the recurrent synaptic weights among HD cells, $W^{\text{HR}}$ is a $N^{\text{HD}}\times N^{\text{HR}}$ matrix of the synaptic weights from HR to HD cells, ${\mathit{r}}^{\text{HR}}$ and ${\mathit{r}}^{\text{HD}}$ are vectors of the firing rates of HR and HD cells respectively, $I_{inh}^{\text{HD}}$ is a constant inhibitory input common to all HD cells, and ${\mathit{n}}^d$ is a random noise input to the axon-distal compartment. ${\mathit{n}}^{\mathit{d}}$ is drawn IID from $N(0,1)$, and its variance is scaled by ${\sigma }_n^2$. Note that in the main text we set ${\sigma }_n$ to zero, but we explore different values for this parameter in Appendix 1. The constant current $I_{inh}^{\text{HD}}$ is in line with a global-inhibition model with local recurrent connectivity, as opposed to having long-range inhibitory recurrent connectivity (Kim et al., 2017). The inhibitory current $I_{inh}^{\text{HD}}$ suppresses HD bumps in general; however the exact strength of this inhibition is not important in our model.

Since several electrophysiological parameters of the fly neurons modeled here are unknown, we use dimensionless conductance values. Therefore, in Equation 2, which describes the dynamics of the axon-distal input of HD cells, currents (e.g. ${\mathit{I}}^d$, $I_{inh}^{\text{HD}}$, and ${\mathit{n}}^d$) are dimensionless. Membrane voltages are also chosen to be dimensionless, and because we measure firing rates in units of 1 /s, all synaptic weights (e.g. $W^{\text{rec}}$ and $W^{\text{HR}}$) then have, strictly speaking, the unit ‘seconds’ (s), even though we mostly suppress this unit in the text. Importantly, all time constants (e.g. ${\tau }_s$), which define the time scale of dynamics, are measured in units of time (in seconds).

Our model incorporates several time scales, whose interplay is not obvious. To facilitate understanding, we summarize the parameters that define the time scales in Appendix 4—table 1, and discuss their relation in Appendix 4.

The axon-distal voltage ${\mathit{V}}^d$ of HD cells is a low-pass filtered version of the input current ${\mathit{I}}^d$, that is,

where ${\tau }_l$ is the leak time constant of the axon-distal compartment. The voltage ${\mathit{V}}^d$ and the current ${\mathit{I}}^d$ have the same unit (both dimensionless), which means that the leak resistance of the axon-distal compartment is also dimensionless, and we assume that it is unity for simplicity. We choose values of ${\tau }_l$ and ${\tau }_s$ (for specific values, see Table 1) so that their sum matches the phenomenological time constant of HD neurons (E-PG in the fly), while ${\tau }_s$ equals to the phenomenological time constant of HR neurons (P-EN1 in the fly, Turner-Evans et al., 2017). Note that ${\mathit{V}}^d$ is the low-frequency component of the axon-distal voltage originating from postsynaptic potentials, that is excluding occasional high-frequency contributions from backpropagating action potentials.

The axon-proximal voltage ${\mathit{V}}^{\mathit{a}}$ of HD cells is then given by

where $C$ is the capacitance of the membrane of the axon-proximal compartment, $g_L$ is the leak conductance, $g_D$ is the conductance of the coupling from axon-distal to axon-proximal dendrites, ${\mathit{I}}^{vis}$ is a vector of visual input currents to the axon-proximal compartment of HD cells, ${\mathit{I}}_{exc}^{\text{HD}}$ is an excitatory input to the axon-proximal compartment, and ${\mathit{n}}^{\mathit{a}}$ is a random noise vector injected to the axon-proximal compartment, drawn IID from $N(0,1)$. The excitatory current ${\mathit{I}}_{exc}^{\text{HD}}$ is assumed to be present only in light conditions. The values of $C$, $g_L$, and $g_D$ in the fly HD (E-PG) neurons are unknown, thus we keep these parameters unitless, and set their values to the ones in Urbanczik and Senn, 2014. Note that since conductances are dimensionless here, $C$ is effectively a time constant.

Following Hahnloser, 2003, the visual input to the i-th HD cell is a localized bump of activity at angular location ${\theta }_i$:

where $M$ scales the bump’s amplitude, $\sigma$ controls the width of the bump, ${\theta }_i$ is the preferred orientation of the i-th HD neuron , ${\theta }_0(t)$ is the position of a visual landmark at time $t$ in head-centered coordinates, and ${\displaystyle I_o^{vis}<0}$ is a constant inhibitory current that acts as the baseline for the visual input. We choose $M$ so that the visual input can induce a weak bump in the network at the beginning of learning, and we choose $\sigma$ so that the resulting bump after learning is ∼60 deg wide. Note that the bump in the mature network has a square shape (Figure 3—figure supplement 2B); therefore we elect to make it slightly narrower than the average full width at half maximum of the experimentally observed bump (∼80 deg; Seelig and Jayaraman, 2015; Kim et al., 2017; Turner-Evans et al., 2017). In addition, the current $I_o^{vis}$ is negative enough to make the visual input purely inhibitory, as reported (Fisher et al., 2019). The visual input is more inhibitory in the surround to suppress activity outside of the HD receptive field. Therefore, the mechanism in which the visual input acts on the HD neurons is disinhibition.

The firing rate of HD cells, which is set by the voltage in the axon-proximal compartment, is given by

where

is a sigmoidal activation function applied element-wise to the vector ${\mathit{V}}^{\mathit{a}}$. The variable $f_{max}$ sets the maximum firing rate of the neuron, $\beta$ is the slope of the activation function, and $x_{1/2}$ is the input level at which half of the maximum firing rate is attained. The value of $f_{max}$ is arbitrary, while $\beta$ is chosen such that the activation function has sufficient dynamic range, and $x_{1/2}$ is chosen such that for small negative inputs the activation function is non-zero.

We note that the saturation of the activation function $f$ in Equation 7 is an essential feature of our model, especially for the convergence of the plasticity rule in Equation 12; see also the section ‘Synaptic Plasticity Rule’. Even though, to the best of our knowledge, it is currently not known whether E-PG neurons actually reach saturation, other Drosophila neurons are known to reach saturation with increasing inputs, instead of some sort of depolarization block (Wilson, 2013; Brandão et al., 2021). Saturation with increasing inputs may be due to, for instance, short-term synaptic depression: beyond a certain frequency of incoming action potentials, the synaptic input current is almost independent of that frequency (Tsodyks and Markram, 1997; Tsodyks et al., 1998).

The firing rates of the HR cells are given by

where ${\mathit{r}}^{\text{HR}}$ is the vector of length $N^{\text{HR}}$ of firing rates of HR cells, the $N^{\text{HR}}\times N^{\text{HD}}$ matrix $W^{\text{HD}}$ encodes the fixed connections from the HD to the HR cells, ${\mathit{r}}_{\text{LP}}^{\text{HD}}$ is a low-pass filtered version of the firing rate of the HD cells where the filter accounts for delays due to synaptic transmission in the incoming synapses from HD cells, ${\mathit{I}}^{vel}$ is the angular velocity input, $I_{inh}^{\text{HR}}$ is a constant inhibitory input common to all HR cells, and ${\mathit{n}}^{\text{HR}}$ is a random noise input to the HR cells drawn IID from $N(0,1)$. We set $I_{inh}^{\text{HR}}$ to a value that still allows sufficient activity in the HR cell bump, even when the animal does not move. The low-pass filtered firing-rate vector ${\mathit{r}}_{\text{LP}}^{\text{HD}}$ is given by

and the angular-velocity input to the i-th HR neuron is given by

where $k$ is the proportionality constant between head angular velocity and velocity input to the network, $v(t)$ is the head angular velocity at time $t$ in units of deg/s, and the factor $q$ is chosen such that the left (right) half of the HR cells are primarily active during leftward (rightward) head rotation. Note that the same ${\tau }_s$ is in both Equation 2 and Equation 9. Finally, as mentioned earlier, the matrix $W^{\text{HD}}$ encodes the hardwired 1-to-1 HD-to-HR connections, i.e., $W_{ij}^{\text{HD}}=w^{\text{HD}}$ if the j-th HD neuron projects to the i-th HR neuron, and $W_{ij}^{\text{HD}}=0$ otherwise. Specifically, for $j$ odd, HD neuron $j$ projects to L-HR neuron $i=\frac{j+1}{2}$, whereas for $j$ even, HD neuron $j$ projects to R-HR neuron $i=30+\frac{j}{2}$. The synaptic strength $w^{\text{HD}}$ is chosen such that the range of the firing rates of the HD cells is mapped to the entire range of firing rates of the HR cells. Specifically, we set $w^{\text{HD}}=\frac{A_{active}}{f_{max}}$, where $A_{active}$ is the range of inputs for which $f$ has not saturated, i.e., the input values for which $f$ remains between about 7% and 93% of its maximum firing rate $f_{max}$ (see Equation 7). Finally, the proportionality constant $k$ is set so that the firing rate of HR neurons does not reach saturation for the range of velocities relevant for the fly (approx. [-500, 500] deg/s), given all other inputs they receive.

In our network, the associative HD neurons receive direct visual input in the axon-proximal compartment and indirect angular velocity input in the axon-distal compartment through the HR-to-HD connections (Figure 1D). We hypothesize that the visual input acts as a supervisory signal that controls the axon-proximal voltage ${\mathit{V}}^{\mathit{a}}$ directly, and the latter initiates spikes. Therefore, the goal of learning is for the axon-distal voltage ${\mathit{V}}^d$ to predict the axon-proximal voltage by changing the synaptic weights $W^{\text{rec}}$ and $W^{\text{HR}}$. This change is achieved by minimizing the difference between the firing rate $f({\mathit{V}}^{\mathit{a}})$ in the presence of visual input and the axon-distal prediction $f({\mathit{V}}^{ss})$ of the firing rate in the absence of visual input. In the latter case and at steady-state, the voltage $V_i^{ss}$ for the i-th HD neuron is an attenuated version of the axon-distal voltage,

with conductance $g_D$ of the coupling from the axon-distal to axon-proximal dendrites and leak conductance $g_L$ of the axon-proximal compartment, as explained in Equation 4, and $p=\frac{g_D}{g_D+g_L}$ in Equation 1. Therefore, following Urbanczik and Senn, 2014, we define the plasticity-induction variable ${\text{PI}}_{ij}$ for the connection between the j-th presynaptic neuron and i-th postsynaptic neuron as

where ${\text{P}}_j$ is the postsynaptic potential of neuron $j$, which is a low-pass filtered version of the presynaptic firing rate r_j. That is,

where * denotes convolution. The transfer function

is derived from the filtering dynamics in Equation 2 and Equation 3 and accounts for the delays introduced by the synaptic time constant ${\tau }_s$ and the leak time constant ${\tau }_l$. In Equation 14, $u(t)$ denotes the Heaviside step function, that is, $u(t)=1$ for ${\displaystyle t>0}$ and $u(t)=0$ otherwise. The plasticity-induction variable is then low-pass filtered to account for slow learning dynamics,

and the final weight change is given by

where $\eta$ is the learning rate and $W_{ij}$ is the connection weight from the j-th presynaptic neuron to the i-th postsynaptic neuron. Note that the synaptic weight $W_{ij}$ is an element of either the matrix $W^{\text{rec}}$ or the matrix $W^{\text{HR}}$ depending on whether the presynaptic neuron $j$ is an HD or an HR neuron, respectively. The value of the plasticity time constant ${\tau }_{\delta }$ is not known, therefore we adopt the value suggested by Urbanczik and Senn, 2014.

Equation 12 is a ‘delta-like’ rule that can be interpreted as an extension of the Hebbian rule; compared to a generic Hebbian rule, we have replaced the postsynaptic firing rate $f(V_i^a)$ by the difference between $f(V_i^a)$ and the predicted firing rate $f(V_i^{ss})$ of the axon-distal compartment of the postsynaptic neuron. This difference drives plasticity in the model. We note that $f(V_i^a)$ is a continuous approximation of the spike train of the postsynaptic neuron, which could be available at the axon-distal compartment via back-propagating action potentials (Larkum, 2013). Furthermore, the axon-distal voltage $V_i^d$ and postsynaptic potentials are by definition available at the synapses arriving at the axon-distal compartment. Note that even though $f(V_i^{ss})$ is the firing rate in the absence of visual input, it can still be computed at the axon-distal compartment when the visual input is available; $V_i^d$ is the local voltage and therefore only a constant multiplicative factor (Equation 11) and the static nonlinearity $f$ need to be computed to obtain $f(V_i^{ss})$. Therefore, the learning rule is biologically plausible because all information is locally available at the synapse.

The learning rule used here differs from the one in the original work of Urbanczik and Senn, 2014 because we utilize a rate-based version instead of the original spike-based version. Even though spike trains can introduce Poisson noise to $f(V_i^a)$, Urbanczik and Senn, 2014 show that once learning has converged, asymmetries in the weights due the spiking noise are on average canceled out.

Another difference in our learning setup is that, unlike in Urbanczik and Senn, 2014, the input to the axon-proximal compartment does not reach zero in equilibrium (see, e.g. Figure 3D, and Appendix 5). Therefore, an activation function with a saturating non-linearity, as in Equation 7, is crucial for convergence, which could not be achieved with a less biologically plausible threshold-linear activation function. This lack of strict convergence in our setup is responsible for the square form of the bump (Figure 3—figure supplement 2B and Appendix 5).

We train the network with synthetically generated angular velocities, simulating head turns of the animal. $W^{\text{rec}}$ and $W^{\text{HR}}$ are both initialized with random connectivity drawn from a normal distribution with mean 0 and standard deviation $1/\sqrt{N^{\text{HD}}}$, as common practise in the modeling literature. In further simulations with various other initial conditions (e.g. in the simulations with gain changes in Figure 4 or in simulations in which we randomly shuffled weights after learning, not shown), we confirmed that the final PI performance is virtually independent of the initial distribution of weights $W^{\text{rec}}$ and $W^{\text{HR}}$.

The network dynamics are updated in discrete time steps $\mathrm{\Delta }t$ using forward Euler integration. The entrained angular velocities cover the range of angular velocities exhibited by the fly, which are at maximum ∼500 deg/s during walking or flying (Geurten et al., 2014; Stowers et al., 2017). The angular velocity $v(t)$ is modeled as an Ornstein-Uhlenbeck process given by

where $\alpha =\mathrm{\Delta }t/{\tau }_v$ and ${\tau }_v$ is the time constant with which $v(t)$ decays to zero, $n(t)$ is noise drawn from a normal distribution with mean 0 and standard deviation 1 at each time step, and ${\sigma }_v$ scales the noise strength.

We pick ${\sigma }_v$ and ${\tau }_v$ so that the resulting angular velocity distribution in Figure 3—figure supplement 2C and its time course, for example in Figure 2A, are similar to what has been reported in flies during walking or flying (Geurten et al., 2014; Stowers et al., 2017). Finally, note that we train the network for angular velocities a little larger than what flies typically display (up to ±720 deg/s).

In Equation 12 we have used the learning error

which controls learning in the i-th associative HD neuron. To quantify the mean learning error $\text{err}(t)$ in the whole network at time $t$, we average $E_i$ across all HD neurons and across a small time interval $[t,t+t_w]$, that is,

with $t_w=10$ s. In Figure 3D, we plot this mean error at every 1% of the simulation, for 12 simulations, and averaged across the ensemble of the simulations. Note that individual simulations occasionally display ‘spikes’ in the error. Large errors occur if the network happens to be driven by very high velocities that the network does not learn very well because they are rare; larger errors also occur for very small velocities, that is, when the velocity input is not strong enough to overcome the local attractor dynamics, as seen, for example, in Figure 2C. On average, though, we can clearly see that the mean learning error decreases with increasing time and settles to a small value (e.g. Figure 3D and Figure 4A).

To decode from the activity of HD neurons an average HD encoded by the network, we use the population vector average (PVA). We thus first convert the tuning direction ${\theta }_i$ of the i-th HD neuron to the corresponding complex number $e^{j{\theta }_i}$ on the unitary circle, where $j$ is the imaginary unit. This complex number is multiplied by the firing rate $r_i^{\text{HD}}$ of the i-th HD neuron, and then averaged across neurons to yield the PVA

The PVA is a vector in the 2-D complex plane and points to the center of mass of activity in the HD network. Finally, we take the angle $\theta$ of the PVA as a measure for the current heading direction represented by the network.

To quantify the variability of heading direction in the trained networks, we define the diffusion coefficient $D$ as:

where $\mathrm{\Delta }\theta$ is the change in heading direction in a time interval $t_{\text{sim}}$. Therefore, $D$ is given by the variance of the distribution of displacements in a given time interval, divided by the time interval.

In the main text, we estimate $D$ during PI, i.e. with velocity inputs only. In this setting, $D$ is the rate at which the variance of the PI errors increases (see e.g. Figure 2B). Deviations from gain-1 PI contribute to this estimate; hence, to single out the effects of noise during training on the stability of the learned attractor in Appendix 1, we also estimate $D$ in the presence of test noise when no inputs are received at all.

Our model assumes the segregation of visual inputs to HD (E-PG) cells from head rotation and recurrent inputs to the same cells. To test this hypothesis, we leverage on the fly hemibrain connectome (Xu, 2020; Clements et al., 2020). First, we randomly choose one E-PG neuron per wedge of the EB, for a total of 16 E-PG neurons. We reasoned this sample would be sufficient because the way E-PG neurons in the same wedge are innervated is expected to be similar. We then find all incoming connections to these neurons from visually responsive ring neurons R2 and R4d (Omoto et al., 2017; Fisher et al., 2019). These are the connections that arrive at the axon-proximal compartment in our model. We then find all incoming connections from P-EN1 cells, which correspond to the HR neurons, and from P-EN2 cells, which are involved in a recurrent excitatory loop from E-PG to P-EG to P-EN2 and back to E-PG (Turner-Evans et al., 2020). These are the connections that arrive at the axon-distal compartment in our model.

To further support the assumption that visual inputs are separated from recurrent and HR-to-HD inputs in the Drosophila EB, we perform binary classification between the two classes (R2 and R4d vs. P-EN1 and P-EN2). We use SVMs with Gaussian kernel, and perform nested 5-fold cross validation, for a total of 30 model runs for every neuron tested (Figure 1—figure supplement 1).

To quantify PI performance of the network and compare to fly performance, we use the measure defined by Seelig and Jayaraman, 2015 and estimate the correlation coefficient between the unwrapped PVA and true heading in darkness. We estimate the correlation in 140 s long trials and report the point estimate and 95% confidence intervals (Student’s t-test, $N=100$).

All code used in this work is available at https://github.com/panvaf/LearnPI (copy archived at swh:1:rev:c6e354f80bf435114e577af70892db41c3ce5315, Vafidis, 2022). The files required to reproduce the figures can be found at https://gin.g-node.org/pavaf/LearnPI.

In the main text, the only source of stochasticity in the network came from the angular velocity noise in the Ornstein-Uhlenbeck process (Materials and methods, Equation 17). Biological HD systems, however, are subject to other forms of biological noise like randomness of ion channels. To address that, we include Gaussian IID synaptic current noise to every location in the network where inputs arrive: the axon-proximal and axon-distal compartments of HD cells and the HR cells (see Materials and methods, parameterized by ${\sigma }_n$ in Equation 2, Equation 4, and Equation 8). We then ask how robustly can the network learn in the presence of such additional stochasticity.

To quantify the network’s robustness to noise, we need to define a comparative measure of useful signals vs. noise in the network. By ‘signals’ we refer to the velocity/visual inputs and any network activity resulting from them, whereas ‘noise’ is the aforementioned Gaussian IID variables. We thus define the signal-to-noise ratio (SNR) to be the squared ratio of the active range $A_{active}$ of the activation function $f$ (defined in Equation 7) over two times the standard deviation of the Gaussian noise, ${\sigma }_n$, i.e.

This definition is motivated by the fact that $A_{active}$ determines the useful range that signals in the network can have. If any of the signals exceed this range, they cannot impact the network in any meaningful way because the neuronal firing rate has saturated, unless they are counterbalanced by other signals reliably present. The factor 2 in the denominator is due to the fact that the noise can extend to both positive and negative values, whereas $A_{active}$ denotes the entire range of useful inputs.

Here, we vary the SNR and observe its impact on learning and network performance. Appendix 1—figure 1A–D shows the performance of a network that has been trained with $\text{SNR}\approx 2$. The resulting network connectivity remains circularly symmetric and maintains the required asymmetry in the HR-to-HD connections for L- and R-HR cells (data not shown). Therefore we plot only the profiles in Appendix 1—figure 1B, which look very similar to the ones in Figure 3C trained with $\text{SNR}=\mathrm{\infty }$. The peak of the local excitatory connectivity in $W^{\text{rec}}$ is not as pronounced. This happens because the noise corrupts auto-correlations of firing during learning.

The network activity still displays a clear bump that smoothly follows the ground truth in the absence of visual input (Appendix 1—figure 1A). There are only minor differences compared to the network without noise (Figure 2A). The presence of the noise is most obvious in the HR cells, since HD cells that do not participate in the bump are deep into inhibition, and therefore synaptic input noise does not affect as much their activity. We note that the network can no longer sustain a bump in darkness when $\text{SNR}=1$, that is when the standard deviation of the noise covers the full active range of inputs (data not shown).

Finally, the neural velocity slightly overestimates the head angular velocity (Appendix 1—figure 1C compared to Figure 2C), and the PI errors diffuse faster in the network with noise ($88.1{\text{deg}}^2/\text{s}$ in Appendix 1—figure 1D compared to $24.5{\text{deg}}^2/\text{s}$ in Figure 2B); these values are also indicated in Appendix 1—figure 1E (triangles). Importantly, we find that the diffusion assumption holds, because the estimation of the diffusion coefficient when varying simulation time (between 10 and 60 s) is consistent.

Appendix 1—figure 1

Download asset Open asset

So far we have addressed the diffusivity of PI in networks that receive velocity input, and we have compared the performance of networks that were trained with and without synaptic input noise. However, in mature networks, that is during testing, it is unclear how large the impact of synaptic input noise is, compared to noise that originates from imperfect PI. To disentangle these two noise contributions during testing, we study diffusivity in networks that do not receive velocity inputs at all (i.e. without PI). In the absence of velocity inputs, we vary the level of synaptic input noise during testing (called ‘test noise’), and estimate diffusion coefficients. Specifically, for each test noise magnitude ${\sigma }_n$ we run 1000 simulations of $t_{\text{sim}}=10$ s, each time randomly initializing the network at one of the angular locations ${\theta }_i$ for which HD neurons are tuned for. For ${\sigma }_n=0$, we run one simulation per ${\theta }_i$, since the simulation is deterministic.

Our results in Appendix 1—figure 1E show that the diffusion coefficients obtained in networks without velocity inputs (dots) are always much smaller that those with velocity inputs (triangles). Thus, imperfect integration of velocity inputs is by far the dominating source of noise in trained, mature networks. Omitting the velocity inputs, we detected small differences between networks that were trained without synaptic input noise (blue dots) and with such noise (orange dots, train noise ${\sigma }_n=0.7$). The network trained with noise is sligthly more diffusive up to the level of test noise at which is was trained (${\sigma }_n=0$ – 0.7), and it is slightly less diffusive beyond that level (${\displaystyle {\sigma }_n>0.7}$).

We conclude that learning PI is robust to synaptic input noise during learning, and that synaptic input noise during testing degrades performance much less than errors due to deviations from perfect gain-1 PI, which are already quite small.

In the main text we trained networks for a set of angular velocities that cover the full range exhibited by the fly (${\displaystyle |v|<500}$ deg/s), and we showed that the mature network can account for several key experimental findings. However, the ability of any continuous attractor network to path-integrate is naturally limited for high angular velocities, due to the synaptic delays inherent in any such network (Zhong et al., 2020). To evaluate the ability of our network to integrate angular velocities, we sought to identify a limit of what velocities could be learned.

The width of the HD bump in our network is here termed BW, and it is largely determined by the width $\sigma$ of the visual receptive field. This is because during training we force the network to produce a bump with a width matching that of the visual input, and this width is then maintained when the latter is not present. The reason for this behavior is that the width of the learned local excitatory connectivity profile in $W^{\text{rec}}$ that guarantees such stable bumps of activity will be similar to the width of the bump, because recurrent connections during learning are only drawn from active neurons (non-zero ${\text{P}}_j$ in Equation 12). As mentioned in the main text, this emphasizes the Hebbian component of our learning rule (fire together — wire together). As a result, the width of local excitatory recurrent connections should be approximately BW.

In Figure 3—figure supplement 1 we show that the higher angular velocities are served by the long-range excitatory connections in $W^{\text{HR}}$. However, these connections might not be strong enough to move the bump by themselves; a contribution from HD cells might still be needed to move the bump at such high angular velocities. In that case, the width of the connectivity bump in $W^{\text{rec}}$ might limit how far away from the current location the bump can be moved. In addition, there is a limitation in how quickly the bump can be moved: the learning rule in Equation 1 tries to predict the next state of the network from the current state; but to activate the next HD neurons in line, current HD and HR cell activity must go through the synaptic delay ${\tau }_s$. Therefore, the maximum velocity that the network can achieve without external guidance (i.e. without visual input) should be inversely proportional to ${\tau }_s$, i.e.

where we assume that $b$ might reflect an effective HD connectivity bump width, which depends on $\sigma$ but also on the angular resolution of the HD network $\mathrm{\Delta }\varphi$, due to discretization effects. In reality, the HR-to-HD connectivity profiles in $W^{\text{HR}}$ likely also have a bearing on $b$.

We then systematically vary ${\tau }_s$ and test what velocities the network can learn. We indeed find that networks can path-integrate all angular velocities up to a limit, but not higher than that. As predicted, this limit is inversely proportional to ${\tau }_s$, for a wide range of delays (Appendix 1—figure 1A). Furthermore, $b$ matches BW reasonably well. Fitting Equation 23 to the data we obtain $b(0.25,6\text{deg})\approx \text{BW}=96\text{deg}$ for $N^{\text{HD}}=N^{\text{HR}}=120$ and $b(0.15,12\text{deg})=75\text{deg}$, $\text{BW}=60\text{deg}$ for $N^{\text{HD}}=N^{\text{HR}}=60$.

As mentioned in the main text, there are two limitations other than synaptic delays why the network could not learn high angular velocities: limited training of these velocities, and saturation of HR cell activity. These limitations kick in for ${\displaystyle {\tau }_s<150\text{ms}}$, for which $v_{max}$ matches the maximum velocity the fly displays (500 deg/s). Therefore to create Appendix 1—figure 1A for these delays, we increased the standard deviation of the velocity noise in the Ornstein-Uhlenbeck process to ${\sigma }_v=800\text{deg/s}$ to address the first limitation, and we increased the dynamic range of angular velocity inputs by decreasing the proportionality constant in Equation 10 to $k=1/540\text{s/deg}$ to address the second.

The velocity gain plot for an example network with high synaptic delays (${\tau }_s=190\text{ms}$) is shown in Appendix 1—figure 1B. Interestingly, we notice that the performance drop at the velocity limit is not gradual; instead, the neural velocity abruptly drops to a near-zero value once past the velocity limit. Further investigation reveals that for velocities higher than this limit, the network can no longer sustain a bump (Appendix 1—figure 1C). This happens because the HD network cannot activate neurons downstream fast enough to keep the bump propagating, and therefore the bump disappears and the velocity gain plot becomes flat.

Equation 23 is similar to a relationship reported in Turner-Evans et al., 2017 (their page 35, 1st paragraph), where it was demonstrated that the phase shift of the HR population bump compared to the HD bump limits angular velocity. Our result hence generalizes this finding in the case where recurrent connections between HD neurons are also allowed.

Finally, we note that so far we only tested the limits of network performance when increasing ${\tau }_s$. To demonstrate that smaller delays also work, as an extreme example we show PI performance in a network where ${\tau }_s=1\text{ms}$ in Appendix 1—figure 1D and E. A potential issue with such small synaptic delays is that the network would not be able to distinguish rightward from leftward rotation for small angular velocities, because the motion direction offset of the HR bumps would be small, and the activity in the two HR populations comparable. In such a setting it is harder to learn the asymmetries in the HR-to-HD connections required to differentiate leftward from rightward movement. Indeed, this effect is visible in Appendix 1—figure 1E where the amplitude of HR-to-HD connections has been suppressed, and in Appendix 1—figure 1D where the flat region for small angular velocities has been extended compared to Figure 2C.

Overall, these results indicate that the network learns to path-integrate angular velocities up to a fundamental limit imposed by the architecture of the HD system in the fly. Furthermore, we conclude that the phenomenological delays observed in the fly HD system in Turner-Evans et al., 2017 are not fundamentally limiting the system’s performance, since they can support PI for angular velocities much higher than the ones normally displayed by the fly.

Appendix 2—figure 1

Download asset Open asset

The networks we have trained in the Results had a circular symmetric initial architecture, including the hardwired HD-to-HR connections $W^{\text{HD}}$. However, such symmetry is unrealistic for any biological system that is assembled by imperfect processes; deviations from symmetry should be expected. Therefore, in this Appendix we let $W^{\text{HD}}$ vary randomly, and observe how PI performance is affected.

First, we remind the reader that the magnitude $w^{\text{HD}}$ of the HD-to-HR connections is chosen so that we take advantage of the full dynamic range of HR neurons; however, the exact magnitude should not be critical for our model. Homeostatic plasticity could adjust the magnitude, but for simplicity we have not incorporated such plasticity rules in our model. Instead, to see whether the exact values of synaptic weights, their circular symmetry, and the 1-to-1 nature of the HD-to-HR connections is crucial for our model, we draw connection strengths randomly.

In a first approach, we let HD neurons project also to adjacent HR neurons. Specifically, if $U(a,b)$ denotes the uniform distribution in the interval $(a,b)$, we sample the magnitude of weights from $U(\frac{w^{\text{HD}}}{2},w^{\text{HD}})$ for the main diagonal and $U(0,\frac{3w^{\text{HD}}}{4})$ for the side diagonals of $W^{\text{HD}}$ (Appendix 3—figure 1A). We then adjust the network connectivity ($W^{\text{rec}}$ and $W^{\text{HR}}$) in a learning phase, similar to the one illustrated in Figure 3E and F. After learning, as shown in Appendix 3—figure 1C-E, PI is still excellent because the learning rule can balance out any deviations from circular symmetry in $W^{\text{HD}}$. It does so by introducing deviations from circular symmetry in the learned weights, mainly in $W^{\text{HR}}$ (Appendix 3—figure 1B). Thus, small deviations from circular symmetry in the learned weights are essential for PI.

Appendix 3—figure 1

Download asset Open asset

We illustrate the necessity to counterbalance small deviations of circular symmetry again in a second example that is based on the symmetric network studied in the Results. Here, we use the connectivity of the network illustrated in Figure 3A–C and also preserve the 1-to-1 nature of HD-to-HR connections, but now we randomly vary their magnitude, while maintaining the same average connection strength; specifically, we sample the magnitude of weights from $U(\frac{w^{\text{HD}}}{2},\frac{3w^{\text{HD}}}{2})$. PI performance in this network is considerably impaired compared to the original network (compare Appendix 3—figure 2 to Figure 2A and C). This is a further argument in favor of synaptic plasticity operating to fine-tune connectivity, because as mentioned we expect that such anatomical asymmetries are indeed present in the biological circuit. Therefore, even if the circular symmetric synaptic weights were passed down genetically with great accuracy, PI performance in flies should be considerably degraded for a biological circuit with anatomical asymmetries when no learning is involved.

Appendix 3—figure 2

Download asset Open asset

To better quantify the effect of anatomical asymmetries, we incorporate both noise in the learned connectivity as in Figure 3—figure supplement 3A,B and noise in the HD-to-HR connections as in Appendix 3—figure 1A. We tune the noise independently for each weight matrix: for $W^{\text{rec}}$ and $W^{\text{HR}}$ we set the variance of the Gaussian noise to $p$ times the variance of the individual weight matrices, while for $W^{\text{HD}}$ we draw the connections connections from $U((1-p)w^{\text{HD}},(1+p)w^{\text{HD}})$ for the main diagonal and $U(0,p{w}^{\text{HD}})$ for the side diagonals. We find that for $p=0.3$ the correlation between the PVA and true heading in darkness drops to $0.27\pm 0.09$ which is below reported fly PI performance (mean correlation across animals $\sim 0.5$ in Seelig and Jayaraman, 2015), while the structure of the weights is preserved. We observed a steep decline in the correlation coefficient between $p=0.25$ (for which the correlation is $0.92\pm 0.04$) and $p=0.3$. Furthermore, we study the effect of perturbing individual weight matrices, and find that perturbing only $W^{\text{HD}}$ with $p=0.3$ considerably affects performance (correlation $0.39\pm 0.09$) while perturbation of $W^{\text{rec}}$ and $W^{\text{HR}}$ together has a much smaller effect on performance. This again argues in favor of learning $W^{\text{rec}}$ and $W^{\text{HR}}$ to counterbalance asymmetries in $W^{\text{HD}}$ (Appendix 3—figure 1). Furthermore, note that confounders other than imperfect weights might be responsible for the degradation of PI performance in the fly, which further argues in favor of learning.

As a final test for the capability of the learning rule to balance anatomical asymmetries, in Appendix 3—figure 3 we use a completely random connectivity for HD-to-HR connections, drawing weights from a folded Gaussian distribution. We find that even then, PI performance of the converged network is great, albeit for a smaller range of velocities. In addition, bumps are not clearly visible in the HR populations anymore; in the main network, HR bumps were inherited from the HD bump due to the sparseness of the HD-to-HR connections. However, when HD-to-HR connections are random, HR cells are no longer mapped to a topographic state space.

Appendix 3—figure 3

Download asset Open asset

We devote this Appendix to discuss requirements for the time scales involved in our model (see Appendix 4—table 1). Several of these time scales are well constrained by biology, and thus we chose to keep them constant. These include the membrane time constants of the axon-proximal and axon-distal compartments, $C/g_L$ and ${\tau }_l$, respectively, which should be in the order of milliseconds; and the velocity decay time constant ${\tau }_v$, for which we choose a value in the same order of magnitude (0.5 s) as experimentally reported (Turner-Evans et al., 2017).

In general, the learning time scale given by $1/\eta$ should be the slowest one in our network model. The time scale should be large enough so that the network samples the input statistics for a long enough time. Varying $1/\eta$ from 2s to 20s to 200s, we find that the final learned weights are virtually identical (not shown). However, $1/\eta$ should not be too large to enable fast enough learning.

The synaptic time constant ${\tau }_s$ is determined from phenomenological delays in the network (Turner-Evans et al., 2017). Nevertheless, we also addressed the impact of varying ${\tau }_s$ in Appendix 2 and found that learning PI was robust in a wide range of values.

In additional simulations, we varied the weight update filtering time constant ${\tau }_{\delta }$ from 0 (effectively removing filtering) to 1000s (which is four orders of magnitude larger than the default value). We observed almost no effect on learning dynamics, and the performance of the final networks was almost identical (not shown). Since the specific value of ${\tau }_{\delta }$ is of little consequence in our network (if $1/\eta$ is large enough), there are hardly any limitations on its value compared to other time scales. Therefore, this justifies ignoring ${\tau }_{\delta }$ in the derivation of a reduced model without noise in Appendix 5.

In this section, we derive a reduced model for the dynamics of the synaptic weights during learning. The goal is to gain an intuitive understanding of the structure obtained in the full model (Figure 3 of the main text). Such a model reduction is obtained by (1) exploiting the circular symmetry in the system; (2) averaging weight changes across different speeds and moving directions; (3) writing dynamical equations in terms of convolutions and cross-correlations. With these methods, we derive a non-linear dynamical system for the weight changes as a function of head direction. Finally, we simulate this dynamical system and inspect how the different variables interact to obtain the final weights. We find that the reduced model results in nearly identical connectivity and learning dynamics to the full network in the main text, and explains how the latter assigns learning errors to the correct weights. Furthermore, it drastically reduces simulation times by two orders of magnitude.

Note that in this section we use slightly different notation compared to the main text. Notably, we refer to the recurrent head direction weight matrix simply as $W$ (omitting the superscript rec), and use capital letters for functions of time and small letters for functions of heading direction.

We study the learning equation (see Equation 12–16 in the main text where the low-pass filtering with time constant ${\tau }_{\delta }$ has been ignored, since we find that the value of ${\tau }_{\delta }$ is not important for learning (see Appendix 4))

where

is the pre-synaptic error at the i-th cell and

is the post-synaptic potential at HD cell $j$, and $H$ is a temporal filter (with time constants ${\tau }_s$ and ${\tau }_l$, see Equation 14 of the main text).

Clockwise movement

Assuming that the head turns clockwise (which equals to rightward rotation, i.e. rotation towards decreasing angles) and anti-clockwise (leftward, i.e. towards increasing angles) with equal probability, we can approximate the weight dynamics by summing the average weight change $W_{ij}^{+}$ for clockwise movement and the average weight change $W_{ij}^{-}$ for anti-clockwise movement:

(27) $\frac{\mathrm{d}}{\mathrm{d}t}{W}_{ij}(t)=\frac{\mathrm{d}}{\mathrm{d}t}{W}_{ij}^{+}(t)+\frac{\mathrm{d}}{\mathrm{d}t}{W}_{ij}^{-}(t).$

We start by assuming head movement at constant speed and we later generalize the results for multiple speeds. We compute the expected weight change $\frac{\mathrm{d}}{\mathrm{d}t}{W}_{ij}^{+}$ for one lap in the clockwise direction at speed ${\displaystyle v^{+}>0}$:

(28) $\frac{\mathrm{d}}{\mathrm{d}t}{W}_{ij}^{+}(t)=\frac{\eta v^{+}}{2\pi }{\int }_0^{2\pi /v^{+}}d\tau E_i^{+}(\tau )P_j^{+}(\tau )$

where

(29) $P_j^{+}(t)={\int }_0^{\mathrm{\infty }}dsH(s)f[V_j^{a+}(t-s)]$

is the post-synaptic potential for clockwise movement, and

(30) $E_i^{+}(t)=f[V_i^{a+}(t)]-f[V_i^{ss+}(t)]$

is the error for a clockwise movement. Assuming that the axon-proximal voltage is at steady state (Equation 4 of the main text with the l.h.s. set to zero and $I_{exc}^{\text{HD}}$ absorbed into $I_{vis}$), the clockwise axon-proximal voltage reads

(31) $V_i^{a+}(t)=V_i^{ss+}(t)+\frac{I_i^{vis}(t)}{g_D+g_L}$

where (see Equation 11 of the main text)

(32) $V_i^{ss+}(t)=\frac{g_D}{g_D+g_L}{V}_i^{d+}.$

From Equation 2 and 3 of the main text, we can write the axon-distal voltage $V_i^{d+}$ as a low-pass filtered version of the total axon-distal current $D_i^{+}$ for clockwise movement (see also Equation 14 of the main text):

(33) $V_i^{d+}={\int }_0^{\mathrm{\infty }}dsH(s)D_i^{+}(t-s),$

which yields

(34) $V_i^{ss+}(t)=\frac{g_D}{g_D+g_L}{\int }_0^{\mathrm{\infty }}dsH(s)D_i^{+}(t-s).$

Importantly, the visual input $I^{vis}$ is translation invariant:

(35) $I_j^{vis}(t)=I_i^{vis}\left(t+\frac{{\theta }_j-{\theta }_i}{v^{+}}\right)$

where ${\theta }_j$ and ${\theta }_i$ are the preferred head directions of the j-th and i-th cell, respectively. As a result of this translation invariance, the recurrent weight matrix $W$ develops circular symmetry:

(36) $W_{ij}=W_{0,(j-i)N_{\text{HD}}}$

where $N_{\text{HD}}$ is the number of HD cells in the system. Consequently, the post-synaptic potential $P_j^{+}$ is also translation invariant:

(37) $P_j^{+}(\tau )=P_i\left(\tau +\frac{{\theta }_j-{\theta }_i}{v^{+}}\right)=P_0(\tau +\frac{\stackrel{{\displaystyle \theta :=}}{\overbrace{{\theta }_j-{\theta }_0}}}{v^{+}}).$

In this case, without loss of generality, we can rewrite Equation 28 for a single row of the matrix $\frac{\mathrm{d}}{\mathrm{d}t}{W}_{ij}^{+}$ as a function of the angle difference $\theta :={\theta }_j-{\theta }_0$:

(38) $\begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}{W}_{ij}^{+}(t)=\frac{\mathrm{d}}{\mathrm{d}t}{W}_{0,(j-i)N_{\text{HD}}}^{+}(t)=\frac{\eta v^{+}}{2\pi }{\int }_0^{2\pi /v^{+}}d\tau E_0^{+}(\tau )P_{(j-i)N_{\text{HD}}}^{+}(\tau )& \end{array}$

(39) $\begin{array}{cc}\hfill =\frac{\eta v^{+}}{2\pi }{\int }_0^{2\pi /v^{+}}d\tau E_0^{+}(\tau )P_0^{+}(\tau +\theta /v^{+})& \end{array}$

(40) $\begin{array}{cc}\hfill =\frac{\eta }{2\pi }{\int }_0^{2\pi }d\phi E_0^{+}(\phi /v^{+})P_0^{+}[(\phi +\theta )/v^{+}]& \end{array}$

(41) $\begin{array}{cc}\hfill =\frac{\eta }{2\pi }{\int }_0^{2\pi }d\phi {ϵ}^{+}(\phi )p^{+}(\phi +\theta )& \end{array}$

(42) $\begin{array}{cc}\hfill =\frac{\eta }{2\pi }({ϵ}^{+}\star p^{+})(\theta )& \end{array}$

(43) $\begin{array}{cc}\hfill =:\frac{\mathrm{d}}{\mathrm{d}t}{w}^{+}(\theta ).& \end{array}$

where we defined ${ϵ}^{+}(\phi ):=E_0^{+}(\phi /v^{+})$ and $p^{+}(\phi ):=P_0^{+}(\phi /v^{+})$, and $\star$ denotes circular cross-correlation.

From Equation 29, we derive

(44) $\begin{array}{cc}\hfill p^{+}(\phi ):=P_0^{+}(\phi /v^{+})={\int }_0^{\mathrm{\infty }}dsH(s)f[V_0^{a+}(\phi /v^{+}-s)]& \end{array}$

(45) $\begin{array}{cc}\hfill \approx {\int }_0^{2\pi }d\beta \underset{=:h^{+}\left(\beta \right)}{\underbrace{\frac{1}{|v^{+}|}H(\beta /v^{+})}}f[\underset{=:v^{a+}\left(\phi -\beta \right)}{\underbrace{V_0^{a+}((\phi -\beta )/v^{+})}}]& \end{array}$

(46) $\begin{array}{cc}\hfill =[h^{+}\ast f(v^{a+})](\phi ).& \end{array}$

The approximation in Equation 45 holds when a bump exists in the network and moves with a velocity below the velocity limit, and it is valid if the temporal filter $H$ is shorter than $2\pi /v^{+}$, that is for $H(t)\ll 1$ for ${\displaystyle t>2\pi /v^{+}}$, which holds for the filtering time constants and velocity distribution we assumed (Appendix 5—figure 1). Therefore, plugging Equation 46 into Equation 43, we obtain:

(47) $\frac{\mathrm{d}}{\mathrm{d}t}{w}^{+}(\theta )\approx \frac{\eta }{2\pi }\{{ϵ}^{+}\star [h^{+}\ast f(v^{a+})]\}(\theta ).$

Appendix 5—figure 1

Download asset Open asset

By using the definition of $ϵ{(\phi )}^{+}$ we derive

(48) ${ϵ}^{+}(\phi ):=E_0^{+}(\phi /v^{+})=f[\underset{=v^{a+}\left(\phi \right)}{\underbrace{V_0^{a+}(\phi /v^{+})}}]-f[\underset{=:v^{ss+}\left(\phi \right)}{\underbrace{V_0^{ss+}(\phi /v^{+})}}]$

with (Equation 31)

(49) $v^{a+}(\phi )=v^{ss+}(\phi )+\underset{=:{\overline{I}}_{vis}\left(\phi \right)}{\underbrace{\frac{I_0^{vis}(\phi /v^{+})}{g_D+g_L}}}$

and (Equation 34)

(50) $\begin{array}{cc}\hfill v^{ss+}(\phi )=\frac{g_D}{g_D+g_L}{\int }_0^{\mathrm{\infty }}dsH(s)D_0^{+}(\phi /v^{+}-s)& \end{array}$

(51) $\begin{array}{cc}\hfill \approx \frac{g_D}{g_D+g_L}{\int }_0^{2\pi }d\beta \underset{=h^{+}\left(\beta \right)}{\underbrace{\frac{1}{|v^{+}|}H(\beta /v^{+})}}\underset{=:d^{+}\left(\phi -\beta \right)}{\underbrace{D_0^{+}\left(\frac{\phi -\beta }{v^{+}}\right)}}& \end{array}$

(52) $\begin{array}{cc}\hfill \approx \frac{g_D}{g_D+g_L}(h^{+}\ast d^{+})(\phi )& \end{array}$

The approximation in Equation 52 is valid if the temporal filter $H$ is shorter than $2\pi /v^{+}$, which again holds true for our parameter choices (Appendix 5—figure 1).

Calculation of the axon-distal input

Let us compute the axon-distal current $D_i^{+}$ to the i-th neuron for clockwise movement. From Equation 2 of the main text, setting the l.h.s. to zero, and splitting the rotation-cell activities in the two populations (L-HR and R-HR), we derive

(53) $D_i^{+}(t)=\underset{:=D_i^{rec+}\left(t\right)}{\underbrace{\sum _j{W}_{ij}(t)f[V_j^{a+}(t)]}}+\underset{:=D_i^{R+}\left(t\right)}{\underbrace{\sum _j{W}_{ij}^R(t)f[V_j^{R+}(t)]}}+\underset{:=D_i^{L+}\left(t\right)}{\underbrace{\sum _j{W}_{ij}^L(t)f[V_j^{L+}(t)]}}+I_{inhib}^{\text{HD}}.$

where $W_{ij}^{\text{R}}$ ($W_{ij}^{\text{L}}$) are the weights from the right (left) rotation cells, and $V_j^{R+}$ ($V_j^{L+}$) are the voltages of the right (left) rotation cells (see Equation 8–10 of the main text):

(54) $\begin{array}{cc}\hfill V_j^{R+}(t)=\frac{A_{\text{active}}}{f_{\text{max}}}{\int }_0^{\mathrm{\infty }}ds{H}_{\text{s}}(s)f[V_j^{a+}(t-s)]+{\overline{I}}_{vel}+I_{inhib}^{\text{HR}}& \end{array}$

(55) $\begin{array}{cc}\hfill V_j^{L+}(t)=& \frac{A_{\text{active}}}{f_{\text{max}}}{\int }_0^{\mathrm{\infty }}ds{H}_{\text{s}}(s)f[V_j^{a+}(t-s)]-{\overline{I}}_{vel}+I_{inhib}^{\text{HR}}.\hfill \end{array}$

The function $H_S(t):=\exp (-t/{\tau }_s)/{\tau }_s$ is a temporal low pass filter with time constant ${\tau }_s$ and the velocity input reads (Equation 10 of the main text)

(56) ${\overline{I}}_{vel}:=v^{+}/(2\pi ).$

Equation 54 and Equation 55 show that the rotation-cell voltages are re-scaled and filtered versions of the corresponding HD-cell firing rates with a baseline shift ${\overline{I}}_{vel}$ that is differentially applied to right and left rotation cells.

From Equation 53, we derive

(57) $d^{+}(\phi ):=D_0^{+}\left(\frac{\phi }{v^{+}}\right)=D_0^{rec+}\left(\frac{\phi }{v^{+}}\right)+D_0^{R+}\left(\frac{\phi }{v^{+}}\right)+D_0^{L+}\left(\frac{\phi }{v^{+}}\right)+I_{inhib}^{\text{HD}}.$

Assuming a large number $N_{\text{HD}}$ of HD cells evenly spaced around the circle, the recurrent axon-distal input reads

(58) $\begin{array}{cc}\hfill D_0^{rec+}\left(\frac{\phi }{v^{+}}\right)={\sum }_j{W}_{0j}f[V_j^{a+}(\phi /v^{+})]+I_{inhib}^{\text{HD}}& \end{array}$

(59) $\begin{array}{cc}\hfill ={\rho }_{\text{HD}}{\int }_0^{2\pi }d\theta w(\theta )f\left[\underset{=:v^{a+}\left(\phi +\theta \right)}{\underbrace{V_0^{a+}\left(\frac{\phi +\theta }{v^{+}}\right)}}\right]+I_{inhib}^{\text{HD}}& \end{array}$

(60) $\begin{array}{cc}\hfill ={\rho }_{\text{HD}}[w\star f(v^{a+})](\phi )+I_{inhib}^{\text{HD}}.& \end{array}$

where ${\rho }_{\text{HD}}=N_{\text{HD}}/2\pi$ is the density of the HD neurons around the circle and we used the fact that the axon-proximal voltage is translation invariant (see also Equation 37):

(61) $V_j^{a+}(\tau )=V_0^{a+}(\tau +\theta /v^{+}).$

Following a similar procedure for $D_0^{R+}$ and $D_0^{L+}$, we obtain:

(62) $d^{+}(\theta )=\left[{\rho }_{\text{HD}}w\star f(v^{a+})+{\rho }_{\text{HR}}{w}^R\star f(v^{R+})+{\rho }_{\text{HR}}{w}^L\star f(v^{L+})\right](\theta )+I_{inhib}^{\text{HD}}$

where ${\rho }_{\text{HR}}=N_{\text{HR}}/2\pi$ is the density of the HR neurons for one particular turning direction (note that we assumed ${\rho }_{\text{HR}}=2{\rho }_{\text{HD}}$ in the main text). In deriving Equation 62 we defined

(63) $\begin{array}{cc}\hfill v^{R+}(\theta ):=V_0^{R+}(t/v^{+})\approx \frac{A_{\text{active}}}{f_{\text{max}}}[h_s^{+}\ast f(v^{a+})](\theta )+{\overline{I}}_{vel}+I_{inhib}^{\text{HR}}& \end{array}$

(64) $\begin{array}{cc}\hfill v^{L+}(\theta ):=V_0^{L+}(t/v^{+})\approx \frac{A_{\text{active}}}{f_{\text{max}}}[h_s^{+}\ast f(v^{a+})](\theta )-{\overline{I}}_{vel}+I_{inhib}^{\text{HR}}.& \end{array}$

where we defined the filter $h_s^{+}(\phi ):=\frac{1}{|v^{+}|}{H}_s(t/v^{+})$, and the approximations are valid if $H_s(t/v^{+})\ll 1$ for ${\displaystyle t>2\pi /v^{+}}$, which holds true for the time constant and velocity distribution assumed in the main text.

Finally, we compute the rotation-cells’ weights change. For these weights, the learning rule is the same as the one for the recurrent connections, except that the post-synaptic HD input is replaced by the post-synaptic HR input. Therefore, following the same procedure as in Equation 38–Equation 46, the rotation weight changes are given by:

(65) $\begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}{w}^{R+}(\theta )=\frac{\eta }{2\pi }\{{ϵ}^{+}\star [h^{+}\ast f(v^{R+})]\}(\theta )& \end{array}$

(66) $\begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}{w}^{L+}(\theta )=\frac{\eta }{2\pi }\{{ϵ}^{+}\star [h^{+}\ast f(v^{L+})]\}(\theta ).& \end{array}$

In summary, for clockwise movement, we obtain the following system of equations:

(67) $\{\begin{array}{cc}{d}^{+}(\theta )\hfill & =[{\rho }_{\text{HD}}w\star f(v^{a+})+{\rho }_{\text{HR}}{w}^R\star f(v^{R+})+{\rho }_{\text{HR}}{w}^L\star f(v^{L+})](\theta )+I_{inhib}^{\text{HD}}\hfill \\ v^{ss+}(\theta )\hfill & =\frac{g_D}{g_D+g_L}(h^{+}\ast d^{+})(\theta )\hfill \\ v^{a+}(\theta )\hfill & =v^{ss+}(\theta )+{\overline{I}}_{vis}(\theta )\hfill \\ v^{R+}(\theta )\hfill & =\frac{A_{\text{active}}}{f_{\text{max}}}[h_s^{+}\ast f(v^{a+})](\theta )+{\overline{I}}_{vel}+I_{inhib}^{\text{HR}}\hfill \\ v^{L+}(\theta )\hfill & =\frac{A_{\text{active}}}{f_{\text{max}}}[h_s^{+}\ast f(v^{a+})](\theta )-{\overline{I}}_{vel}+I_{inhib}^{\text{HR}}\hfill \\ {ϵ}^{+}(\theta )\hfill & =f[v^{a+}(\theta )]-f[v^{ss+}(\theta )]\hfill \\ \frac{\mathrm{d}}{\mathrm{d}t}{w}^{+}(\theta )\hfill & =\frac{\eta }{2\pi }\{{ϵ}^{+}\star [\underset{=:p^{+}}{\underbrace{h^{+}\ast f(v^{a+})}}]\}(\theta )\hfill \\ \frac{\mathrm{d}}{\mathrm{d}t}{w}^{R+}(\theta )\hfill & =\frac{\eta }{2\pi }\{{ϵ}^{+}\star [\underset{=:p^{R^{+}}}{\underbrace{h^{+}\ast f(v^{R+})}}]\}(\theta )\hfill \\ \frac{\mathrm{d}}{\mathrm{d}t}{w}^{L+}(\theta )\hfill & =\frac{\eta }{2\pi }\{{ϵ}^{+}\star [\underset{=:p^{L+}}{\underbrace{h^{+}\ast f(v^{L+})}}]\}(\theta ).\hfill \end{array}$