Ring Attractor Dynamics Emerge from a Spiking Model of the Entire Protocerebral Bridge

Animal navigation is accomplished by a combination of landmark-following and dead reckoning based on estimates of self motion. Both of these approaches require the encoding of heading information, which can be represented as an allocentric or egocentric azimuthal angle. Recently, Ca2+ correlates of landmark position and heading direction, in egocentric coordinates, were observed in the ellipsoid body (EB), a ring-shaped processing unit in the fly central complex (CX; Seelig and Jayaraman, 2015). These correlates displayed key dynamics of so-called ring attractors, namely: (1) responsiveness to the position of external stimuli; (2) persistence in the absence of external stimuli; (3) locking onto a single external stimulus when presented with two competitors; (4) stochastically switching between competitors with low probability; and (5) sliding or jumping between positions when an external stimulus moves. We hypothesized that ring attractor-like activity in the EB arises from reciprocal neuronal connections to a related structure, the protocerebral bridge (PB). Using recent light-microscopy resolution catalogs of neuronal cell types in the PB (Lin et al., 2013; Wolff et al., 2015), we determined a connectivity matrix for the PB-EB circuit. When activity in this network was simulated using a leaky-integrate-and-fire model, we observed patterns of activity that closely resemble the reported Ca2+ phenomena. All qualitative ring attractor behaviors were recapitulated in our model, allowing us to predict failure modes of the putative PB-EB ring attractor and the circuit dynamics phenotypes of thermogenetic or optogenetic manipulations. Ring attractor dynamics emerged under a wide variety of parameter configurations, even including non-spiking leaky-integrator implementations. This suggests that the ring-attractor computation is a robust output of this circuit, apparently arising from its high-level network properties (topological configuration, local excitation and long-range inhibition) rather than fine-scale biological detail.


INTRODUCTION
An animal navigating in its environment relies on landmarks to estimate its orientation and position (Collett and Graham, 2004). However, in the absence of visual cues, many animals maintain a representation of their heading and position without landmarks by continuously tracking their own motion to calculate navigation vectors to return to a specific location, a process called path integration (Etienne and Jeffery, 2004). Numerous studies have identified patterns of neural activity that could represent heading, one of the elements needed for path integration. These studies have further shown that heading representations are tuned by visual information but can be updated in the dark, without any visual feedback (Taube, 2007;Seelig and Jayaraman, 2015;Varga and Ritzmann, 2016), presumably by exploiting self-generated motion cues like efference copy (Kim et al., 2015). By integrating heading and distance traveled, an animal can estimate its current position (McNaughton et al., 2007) and calculate a return vector. Heading estimation requires the tracking of variables in angular coordinates, a computation that can be accomplished by ''ring attractor networks'' (Skaggs et al., 1995;Zhang, 1996;Solovyeva et al., 2016).
In theoretical models of ring attractor networks, neighboring nodes connect to form a topological ring. The value of an angular variable is encoded in the radial position of a ''bump'' of neural activity within this ring. This bump arises through the combined dynamics of short range excitation and global or long range inhibition between nodes (Skaggs et al., 1995;Zhang, 1996;Knierim and Zhang, 2012). Asymmetric excitation of adjacent nodes causes the bump to move in the direction of the excitation as its previous position is inhibited. Importantly, these neighboring nodes do not have to be physically adjacent but adjacent only in connectivity. In mammals, ring attractors are thought to explain the dynamics of the head direction (HD) cells, which are primarily found in the thalamus and cortical areas associated with the hippocampus (Taube, 2007). Each HD cell is tuned to a particular head orientation in the horizontal plane and the direction in which the cell fires maximally is referred to as its preferred direction. Activity in topologically neighboring HD cells encodes continuously varying allocentric directions and the motion of the bump through the functionally (though not physically) ring-shaped network of HD cells encodes head orientation. Studying dynamics in mammalian neurons encoding HD is difficult as they are spread across relatively large areas of the brain and not spatially organized according to their preferred directions, making simultaneous monitoring of their activity challenging.
In insects, it was recently shown that a physically ring-shaped concentration of neuronal connections (neuropil) may function as a ring attractor (Seelig and Jayaraman, 2015) within the midline-spanning central complex (CX), of Drosophila melanogaster. Specifically, the ellipsoid body (EB) and protocerebral bridge (PB) may contain a neural circuit implementing a ring attractor. The EB neuropil has a closed ring shape in dipteran insects but is split ventrally and therefore roughly linear or bean-shaped in all other insect groups, where it is called the lower division of the central body (Strausfeld, 1976). Due to the evolutionary conservation of morphological cell types in the CX, it likely retains ring-shaped functional connections in all insects (Pfeiffer and Homberg, 2014). Furthermore, the lower division of the central body has been shown to encode the angular position of the sun in locusts, a continuous variable in angular coordinates, suggesting ring-like function without closed ring shape (Heinze and Homberg, 2007;Homberg et al., 2011;Heinze, 2014). In dipterans, the compact size and physical ring shape of this neuropil uniquely facilitates the study of putative ring attractor dynamics in a complete and intact circuit that can be simultaneously imaged in an awake behaving fly (Seelig and Jayaraman, 2015).
In a closed-loop behavioral setup, Ca 2+ activity in putative dendritic processes of one neuronal population within the EB was shown to encode relative angular position of a vertical stripe on a 2-D light emitting diode screen (Seelig and Jayaraman, 2015). Seelig and Jayaraman (2015) noted several features of their circuit that are typical of ring attractor networks (Haferlach et al., 2007;Knierim and Zhang, 2012;Arena et al., 2013). The Ca 2+ activity in so-called ''wedge neurons'' was localized in a single bump at any one time and this bump moved in response to the animal turning on a ball, mimicking a change in heading. Furthermore, the bump exhibited spatial stability. When the fly was stationary for a long period of time, the bump would sometimes fade, but it typically reappeared in the same location when the fly resumed moving. This suggests the bump can be stored in ways other than intracellular Ca 2+ concentration, such as subthreshold voltages. The bump locked onto a single stripe when two competitor stripes were presented and was observed to jump between identical stripes from time to time. While these bump dynamics were recorded in the context of a closed-loop behavioral assay, causal relationships between the bump, behavior, motor commands and efference copy signals have not yet been identified.
The neurons exhibiting these ring attractor-like dynamics connect two of the neuropil that make up the CX, tiling the EB with dendritic arbors and the PB with presynaptic boutons (Figure 1). They are called E-PGs (EB-PB-Gall neurons, called PBG1-8.b-EBw.s-D/Vgall.b in Wolff et al. (2015), EB.w.s or ''wedge neurons'' in Seelig and Jayaraman (2015) and EIP in Lin et al. (2013); see Table 1 for all abbreviations), denoting the flow of information within them from the EB to the PB and Gall (a secondary structure immediately outside the CX). The EB and PB are notable for their division into columnar segments, known as glomeruli in the PB and wedges/tiles (Wolff et al., 2015) in the EB. These neuropil contain many different neural cell types beyond those shown by Seelig and Jayaraman (2015) to encode angular position. In the PB, these have been recently characterized at the level of morphology using single-cell stochastic labeling methods (Lin et al., 2013;Wolff et al., 2015). The resulting catalog revealed that of the approximately 18 classes of neurons within the PB, only three reciprocally connected the EB and the PB.
We sought to test the hypothesis that PB neurons have circuit dynamics consistent with a ring attractor, using their connectivity as enumerated in these recent mapping articles. With a leaky integrate-and-fire model and simple connectivity rules, derived from light-microscopy resolution neuronal morphologies, we have found that a simple model recapitulates the bump of Ca 2+ activity and essentially all of the in vivo dynamics previously observed (Seelig and Jayaraman, 2015). Furthermore, we have found that this circuit is robust to variation in synaptic weights, behaving like a ring attractor under a wide variety of parameters, perhaps indicating that computing FIGURE 1 | The protocerebral bridge (PB) neural circuit. (A) Diagram of the PB and ellipsoid body (EB), illustrating three out of four modeled neural subtypes, the E-PGs, P-ENs and Pintrs. Not shown are the P-EGs which project from the PB to the EB. Axonal arbors are indicated with circular varicosities/boutons. Dendritic arbors are intricate with fine linear branches. Overlap of an axonal arbor and a dendritic arbor within a single anatomical compartment (gray regions) is sufficient to postulate a synapse between neurons. Neurons with identical morphologies at the level of these anatomical compartments (e.g., the two dark blue E-PGs) are represented in the model as a single neuron. (B) Matrix representation of the connectivity of the PB-EB circuit. A filled rectangle in row i, column j indicates a synapse, with neuron i presynaptic and neuron j postsynaptic. Different fill colors indicate different synapse classes, whose within-class strengths are shown at right. (C) Graph, with node positions determined by a force-directed algorithm of the network with connectivity shown in (B), which forms a ring with bilateral symmetry. Thick edges indicate lateral and reciprocal excitatory loops (local excitation) from neuron 38 (as an example) as well as excitatory connections to inhibitory neurons that target all glomeruli (long-range inhibition). (D) Schematic of key circuit motifs in the PB-EB network. Green and blue arrows represent excitatory neurons, red inhibitory neurons. All E-PGs are shown to illustrate the mapping of the EB to each PB hemisphere. Only one P-EG and one P-EN are shown. The pool of Pintrs are represented as a single neuron. Stars indicate the position where the bump will emerge in the EB and the PB hemispheres. Numbers on labels indicate the sequence in which neurons are activated.

MATERIALS AND METHODS
Simulations were run in MATLAB 2015a and 2016a (The MathWorks, Inc., Natick, MA, USA) using custom scripts. All code to recapitulate these results is available at: http://lab.debivort.org/protocerebral-bridge-ring-attractormodel.

Network Construction
To construct a circuit model of the PB, we began with the catalog of morphologically defined cell types in the PB (Wolff et al., 2015). This work enumerates all neuronal cell types within the PB, characterizing two cells as belonging to the same type if their presynaptic and postsynaptic arbors (as determined by MultiColor FlipOut imaging (Nern et al., 2015)) are in the same neuropil compartments. Compartments are defined as spatially distinct regions of the major glia-ensheathed neuropils of the CX and associated regions. For example, the PB itself contains 18 glomerular compartments and the EB contains 16 wedge and 8 tile compartments. We included in our model: (1) any neuron with postsynaptic processes in the PB and presynaptic processes in other compartments (output neurons), provided there is a PB input neuron with a postsynaptic arbor overlapping the presynaptic arbor of that output neuron; and (2) any neuron with presynaptic processes in the PB and postsynaptic arbors that overlap presynaptic arbors of neurons projecting out of the PB (input neurons; Figure 1). This includes all the neuronal cell types cataloged in Wolff et al. (2015) except for five classes of fan-shaped body (FB) projecting neurons (output only) and two classes of PB input neurons from the posterior slope (input only). We assumed that all neurons could be cleanly divided into dendritic and axonal compartments, and that information flows exclusively from the former to the latter.
The broad classes of neurons that met this criterion were the P-ENs (PB output neurons with axons in the EB and noduli), P-EGs (PB output neurons projecting to the EB and Gall), E-PGs (PB input neurons with dendrites in the EB and output to the Gall), and Pintrs (PB intrinsic neurons with both dendritic arbors and presynaptic boutons in the PB). See Table 1 for abbreviations. P-ENs and P-EGs comprise 16 types each, defined by which PB glomerulus contains their dendrites. E-PGs comprise 18 types, defined by which PB glomerulus contains their axons (unlike ''wedge'' neurons (Seelig and Jayaraman, 2015) E-PGs also include neurons innervating the first and last glomeruli of the PB G9L and G9R (Wolff et al., 2015)). The Pintrs comprise 10 types, defined by which PB glomerulus contains their axons. If their projections were identical at the level of the 60 types described above, individual neurons were considered identical, and represented by a single neuron in the model. Lastly, we assumed that neurons formed no autapses. The connectivity of the network thus defined is shown in Figure 1B. We examined the topological arrangement of this network by using a force-directed algorithm (Fruchterman and Reingold, 1990; Webb and Stone, personal communication) to arrange nodes representing neurons. The connectivity present in this network has a ring-like topology, with bilateral symmetry ( Figure 1C). From this representation, key circuit motifs can be discerned ( Figure 1D). Input depolarizing the E-PGs has the potential to: (1) activate a excitatory E-PG/P-EG loop; (2) activate adjacent E-PGs via the P-ENs; and (3) broadly inhibit all PB glomeruli via the pool of Pintrs.

Circuit Physiological Assumptions
Circuit dynamics were implemented using leaky-integrate-andfire (Stein, 1967) neuronal models, with values for the membrane capacitance, resistance, resting potential, undershoot potential, and postsynaptic current (PSCs) time constants and magnitudes, chosen to reflect generic neuronal properties (Hodgkin and Huxley, 1952). The important free parameters of the model were the strengths and signs of the synapses between each type of neuron. We assume that the strengths of all synapses between two classes of neurons (a ''synapse class,'' e.g., all synapses between P-ENs and E-PGs) were identical.
Strength of synapses was implemented as the number of PSC equivalents per action potential. Excitatory neurons induced positive, depolarizing currents in their postsynaptic partners and inhibitory neurons negative currents. We assumed all neurons were excitatory unless we had evidence otherwise. The Pintrs are glutamatergic (Daniels et al., 2008) and possess connectivity similar to other inhibitory local neurons in spatially compartmentalized neuropils, e.g., the antennal lobe (Chou et al., 2010) and lateral horn (Fis¸ek and Wilson, 2014), therefore we assumed they are inhibitory (Liu and Wilson, 2013).
To deliver inputs to the circuit, we assumed that information first flows into the EB (this assumption has no bearing on our qualitative conclusions). Therefore, for each run of the model, the timing of action potentials in not-explicitly-simulated neurons upstream of the E-PGs was determined. These action potentials induced in the E-PGs excitatory currents with a strength equivalent of one PSC each. We assumed that background activity in these upstream neurons produced a Poisson-process sequence of action potentials with a mean rate of 5 Hz. On top of this, Poisson-process spikes at higher rates (peaking at 120 Hz) in subsets of E-PG types represented sensory-like input into the PB (Figure 2A), e.g., the azimuthal angle of light polarization (Heinze, 2014;Bockhorst and Homberg, 2015) or the retinotopic position of a landmark Jayaraman, 2013, 2015).

Circuit Simulations
The circuit network structure was coded from the morphology descriptions in Wolff et al. (2015) per the rules described in the ''Results'' Section below. Leaky-integrate-and-fire dynamics were used to simulate spiking neurons, and implemented using Euler's method to evaluate the following equation, with ∆t = 10 −4 s: where V i is the membrane voltage of neuron i, I in is input current from neurons outside the PB-EB circuit (0 in all neurons other than the E-PGs), M j,i is the network connectivity matrix with entries equal to the synapse strength (in units of excitatory or inhibitory PSCs), I j is the output current of other neurons in the PB-EB circuit, and I ect is simulated ectopic current (such as might be induced by thermogenetic or optogenetic manipulation). We used parameter values that correspond to a generic spiking neuron (Hodgkin and Huxley, 1952), but these values are consistent with various Drosophila measurements or measurements of PB neurons in other species. C m is the membrane capacitance (0.002 µF in all neurons, assuming a surface area of 10 −3 cm 2 ; Gouwens and Wilson, 2009), V 0 is the resting potential (−52 mV in all neurons; Rohrbough and Broadie, 2002;Sheeba et al., 2008), R m is the membrane resistance (10 M in all neurons; Gouwens and Wilson, 2009). When a neuron's voltage reached the firing threshold of −45 mV (V thr ; Sheeba et al., 2008;Gouwens and Wilson, 2009), a templated action potential trace was inserted into its voltage time series. This trace was defined as follows: Where V max is the (purely cosmetic, as it does not affect the circuit dynamics) peak action potential voltage (20 mV; Rohrbough and Broadie, 2002), V min is the spike undershoot voltage (−72 mV; Nagel et al., 2015), t AP is the length of an action potential (2 ms; Gouwens and Wilson, 2009;Gaudry et al., 2013), Normpdf (a, b, c) is the probability density function of a Gaussian with mean a and standard deviation b at c, and α, β, γ and δ are normalization parameters so that the max and min of the Normpdf and sin segments are 1 and 0 respectively prior to scaling by the voltage terms. The firing of an action potential also triggered the addition of a templated PSC trace to the output current time series of the firing neuron. The PSC trace was defined as follows in terms of t in ms: Where I PSC is the amplitude of a PSC (5 nA; Gaudry et al., 2013); excitatory and inhibitory PSCs were assumed to have the same magnitude but opposite sign), t PSC is the half-life of PSC decay (5 ms; Gaudry et al., 2013), and α , β , γ and δ are normalization parameters so that the max and min of the sin and exponential terms are 1 and 0 respectively prior to scaling by I PSC . PSC traces had a length equal to 2 + 7t PSC ms, corresponding to 2 s of rise time plus 7 times the decay half-life.
Synapse strength parameters were explored manually to identify the baseline configuration in Figure 1. Thereafter, parameter exploration was conducted as described in the ''Results'' Section. The overall magnitude of the synapse strength parameters shown in Figure 1 was the main free parameter of the model. The average synapse strengths of each synapse class are also free parameters, though we found that adjusting only the strengths of the Pintr > P-EG and Pintr > P-EN synapse classes was sufficient to recapitulate bump dynamics.
Leaky-integrator dynamics were used to simulate non-spiking graded potential neurons, and were implemented using Euler's method to evaluate the following equation, with ∆t = 10 −4 s: Where all variables and constants are as defined above, and I max is the maximum PSC achievable in a synapse of strength one within the PB-EB circuit. First, the scaling parameter of the current-voltage tanh transfer function (20) was determined empirically. This value yielded dynamics that were the most bump-like, given the synapse strength parameters determined in the leaky-integrate-and-fire model. Then, synapse strength parameters that produced a maximally functional bump were identified by adding Gaussian noise to the baseline synapse strength parameters from the leaky-integrate-and-fire model. This noise had mean of zero and a standard deviation of 100% of the baseline value of each synapse parameter. The dynamics of approximately 200 such random configurations were examined manually, and those producing the best bump-like behavior were then iteratively refined using them as a new baseline, and then adding Gaussian noise with a standard deviation of 10% and then 5% of the respective baseline values. The formula for these dithered synapse strength parameters was (s k is the strength of synapse class k, s' k is that parameter after dithering): In order to break initial symmetry and allow the bump to move ''spontaneously'' random Gaussian noise with mean zero and standard deviation of 3 × 10 −10 V was added to each neuron in each timestep.
Bump position was estimated and visualized by convolving the action potential rasters of each neuron with a Gaussian kernel with a standard deviation of 24 ms. This approximates a Ca 2+ signal in these neurons. Bump position was determined by taking the centroid, modulo eight, of this convolved representation for the P-ENs in each hemisphere. The estimated centroid of each hemisphere's P-ENs was averaged to produce the final centroid estimate.

Circuit Dynamics Parameter Sensitivity Analysis
Circuit dynamics were captured for multidimensional analysis by simulating the network for 2 s, with inputs representing a rotating bar and two static competitors (setting parameters SweepBarBool and TwinBarBool equal to one in PBexperiment.m). Starting from the baseline synapse strength parameters, we added Gaussian noise to each parameter using the dithering formula above. Circuit dynamics during the rotating bar and competitor bar stimuli are particularly diagnostic of bump performance. We trimmed the circuit ''Ca 2+ signal'' time series to 200 ms from each of these two stimulus periods, then subsampled the time points 20:1, yielding a final 200 diagnostic time points. Dynamics from 10,000 networks with randomly dithered synapse strength parameters were computed, each yielding a single point in a 12,000 dimensional space representing circuit dynamics (200 diagnostic time points × 60 neurons). In addition, we added points representing the dynamics of networks in which a single synapse class parameter value was swept systematically from −9x to 10x its original value. The dynamics from these systematic sweeps were added to the dynamics from the randomly dithered networks and projected into two dimensions using principle components analysis (PCA) for visualization. Clusters of dynamics were enumerated using k-means clustering in the original 12,000 dimensional space. Representative dynamics of each cluster were computed by averaging all of the Gaussian-convolved spike rasters receiving each k-means cluster label.

Bump Circuit Dynamics
As a starting point for our characterization of circuit dynamics, we assumed, rather arbitrarily, that all synapse classes had a strength of 20 PSCs. With a small amount of manual parameter searching, we found that if the inhibitory synapses between the Pintrs and the P-ENs and those between the Pintrs and P-EGs had strengths of 15, circuit activity recapitulated several key phenomena that have been observed in Ca 2+ recordings of the E-PGs (Figures 2B,C; Movie 1): (1) a stable ''bump'' of activity appeared at one position in the glomerular axis of the PB and the corresponding EB position, as observed by Seelig and Jayaraman (2015). This bump was almost always distributed over two or three glomeruli (25%-38% of the azimuthal axis), corresponding roughly to the size of the Ca 2+ bump they imaged (Seelig and Jayaraman, 2015). (2) The bump jumped or slid to the position of a novel sensory cue (i.e., a vertical bar), represented as an increased firing rate in the neurons upstream of a single E-PG. (3) When the position of this input activity moved across adjacent glomeruli (moved in its azimuthal position), the bump followed. (4) When two competing vertical bars were provided in the form of firing-rate-matched activity upstream of two non-adjacent E-PGs, the bump moved to the position of one of the cues. (5) Occasionally, during the presentation of competitor bars, the bump would switch positions from one cue to the other. These characteristics were present for a wide range of synapse strength parameters (see stability analysis below).
As reported by Seelig and Jayaraman (2015), the bump appeared to be fairly stable in the dark (i.e., with only baseline background activity present upstream of the input neurons).
Our baseline synapse strength parameter values yielded a bump ''spontaneous drift rate'' comparable to those observed in vivo (approximately 1 glomeruli/s; Figure 2D). We observed that the angular position encoded by the position of the bump had a highly discretized distribution while drifting in the dark ( Figure 2E); the vast majority of the time, the bump was present in one of 15 azimuthal positions, and among these, ±5π/8 was the most abundant, followed by ±π/8. The distribution of bump speed during spontaneous motion (i.e., any motion in the dark) was trimodal ( Figure 2F). These modes may correspond to staying in position, sliding between adjacent positions, and jumping between non-adjacent positions.

Bump Relation to Action Potentials
The emergence of a bump was remarkably robust, even a single action potential upstream of a single input E-PG was sufficient to induce a bump at the position of that action potential that would persist indefinitely in the absence of further action potentials (Figures 3A,B; Movie 1). We observed that occasionally the bump, as encoded by action potentials, would disappear briefly (up to a few tens of milliseconds at a time; Figures 3C-E). During these periods, none of the PB neurons would fire any action potentials, even if there were occasional action potentials in the neurons upstream of the E-PGs. This implies that the bump can be ''stored'' in sub-threshold potentials. These brief disappearances tended to happen when the bump was located at one of the less frequent azimuthal positions (e.g., ±7π/8).
Several sets of neurons appeared to fire synchronously in the circuit (Figure 3E), specifically, those Pintrs that have axonal arbors in two PB glomeruli, bilaterally paired P-ENs and P-EGs, and bilaterally paired E-PGs (though this group of neurons is somewhat less synchronous by virtue of their being the input neurons that are stimulated at random times by upstream neurons). Leaky-integrator implementations (without action potentials) of this model could also produce a bump that persisted in the absence of sensory input,  Figure 2B). (C) Circuit input via the E-PGs corresponding to 250 ms of darkness. (D) Corresponding dynamics of all neurons in the circuit revealing a ∼20 ms window (shaded gray) in which the bump disappears, i.e., is not represented in action potentials, but reappears in the same position after the window (axes and color scale as in Figure 2B). (E) Corresponding voltage traces. For clarity, the trace of every other neuron has been removed. (F) Depolarizing currents representing input into the E-PGs in a leaky-integrator implementation of the circuit, vs. time. Synapse strength parameters used were those that provided dynamics most closely approximating a bump. (G) Corresponding voltages in the entirety of the circuit. Symbols indicate elements of canonical bump phenomenology, as in Figure 2C. * Indicates "spontaneous" shifts in the position of the bump in darkness. † Indicates the bump jumping to the position of a bar as soon as it appears, and then following it as it rotates. ‡ Indicates the bump jumping to the favoring one of two competitors (the lower competitor, most clearly discernible in the left hemisphere-the top half of each neuron type). selected between competitor bars or formed a unitary bump after competitors were removed (Figures 3F,G). However, the bump in this implementation did not have the same rapid spontaneous bump formation, spatial precision, or strong selectivity between competitors seen in the leaky-integrate-andfire implementation (though it did have weak selectivity between competitors).

Synapse Strength Parameter Sensitivity
We next examined whether we had serendipitously found the only synapse strength parameters that recapitulated so many experimental bump phenomena or were these dynamics robust to parameter values. We added random, Gaussian-distributed noise (mean = 0, standard deviation = 20% of each parameter's baseline value) to the synapse strength parameters and then stimulated these dithered circuits with inputs of: (1) sequential bursts of activity in adjacent wedges representing a rotating bar; and (2) elevated activity in two non-adjacent glomeruli representing stationary competitor bars ( Figure 4A). For each of these configurations, the ensuing circuit activity in all neurons during diagnostic periods of this stimulation (200 ms from the rotating bar phase and 200 ms from the beginning of the competitor bars) were treated as points in a high dimensional space of circuit behavior. These points were clustered and averaged within a cluster to provide an exhaustive catalog of the modes of dynamics that this circuit topology can produce (Figures 4A,B). Of 15 modes, three feature sets of neurons with essentially no activity (modes 1-3) and five feature sets of seizing neurons (modes 11-15). Two of the remaining modes feature bumps exhibiting all the key properties observed experimentally (i.e., those shown in Figure 2): modes 4 and 5, which are distinguished largely by which competitor they select. Two modes have bumps that are stable on too-long timescales and extend over too many adjacent glomeruli, but otherwise show Distributions of synapse parameters that support proper bump function. Gaussian noise from the solid gray distribution was sampled and added to each synapse class independently. This was repeated 24,000 times and the resulting circuit dynamics were k-means clustered into 400 clusters. Six clusters were identified that had proper bump dynamics comparable to modes 4 and 5 of (A). The distribution of synapse strength parameter offsets represented in these clusters is shown for each synapse class strength parameter.
the key properties (modes 6 and 7). The remaining three modes (8-10) have some key bump properties, but are stable on too-long timescales, are too wide, and fail to select between competitor bars. Thus, it appears that bumps with the properties observed by Seelig and Jayaraman (2015) are a robust output of circuits with this topology under a wide range of synapse strength parameters.
To understand the contribution of each synapse class to circuit function, we systematically varied the strength of each synapse class from −9x to 10x its original value (Figure 4C). Converting excitatory drive from the PB to the EB into inhibition (by reversing the sign of either the P-EN > E-PG or P-EG > E-PG synapses) eliminated input-independent bump activity in the P-EGs (mode 3). Increasing the strength of that excitatory drive led to too-stable bumps without competitor selectivity (modes 9 and 10) and eventually seizure across the circuit (mode 11). Increasing inhibition of P-ENs (by either reversing the excitatory E-PG > P-EN synapses or amplifying the strength of the inhibitory Pintr > P-EN synapses), not surprisingly, eliminated activity in the P-ENs (mode 1). Conversely, the opposite manipulations resulted in a too-stable bump (mode 10) and eventually seizure of the P-ENs (and E-PGs and Pintrs; modes 12 and 15). Increasing inhibition of P-EGs (by either reversing the excitatory E-PG > P-EG synapses or amplifying the strength of the inhibitory Pintr > P-EG synapses), not surprisingly, eliminated input-independent bump activity in the P-EGs (mode 3). Conversely, the opposite manipulations resulted in a too-stable bump (mode 8) and eventually seizure of the P-EGs (and E-PGs and Pintrs; mode 13). Increasing inhibition of Pintrs (by either reversing the excitatory E-PG > Pintr synapses or amplifying the strength of the inhibitory Pintr > Pintr), resulted in too-stable bumps (mode 7), bumps with no competitor selectivity (mode 10) and eventually seizure in P-ENs, P-EGs, and E-PGs (mode 14). The opposite manipulations eliminated input-independent bump activity in the P-ENs (mode 2) and eventually all activity in P-ENs and P-EGs (mode 1).
This systematic variation of synapse class strength parameter values also provides evidence of the robustness of the bump phenomenon in this circuit. Increasing or decreasing the strength of a synapse class by up to 50% of its baseline value, for example, seldom changes the mode of circuit dynamics ( Figure 4C). Thus, it seems a sizable parameter subspace around the baseline values can produce bump phenomena. This analysis allows us to assess how much of the parameter space around the baseline produces bumps, not how much of the total space can produce bumps. To discover more distant parameter configurations that might also work, we added a substantial amount of Gaussian noise to all parameters simultaneously (mean = 0, standard deviation = 100% of each parameter's baseline value; this reverses the sign of a parameter 16% of the time). The vast majority of circuits with these more broadly-sampled random parameter configurations seized or were silent in at least one neural subtype but a small portion (∼1.5% of 24,000 random parameter configurations) exhibited correct bump phenomena. These were identified and pooled by k-means clustering of circuit dynamics.
The distributions of each parameter within this pool of bump producing circuits are shown in Figure 4D. Each parameter can evidently take on a wide range of values, and with the right corresponding changes in other parameters, support bump function. Notably, almost all parameters could even have their sign reversed from excitatory to inhibitory or vice versa, and still contribute to a bump-producing circuit. The exceptions were the Pintr > P-EN and Pintr > P-EG synapses, which could be silenced but not converted into excitatory synapses, and still produce a bump. In general, however, the random noise that was added to the parameters in bump-producing circuits was positive, meaning that excitatory synapses could generally be made more excitatory, and inhibitory synapses more inhibitory, without loss of bump function. Several parameter distributions appeared to be multimodal, notably P-EG > E-PG, E-PG > P-EG, E-PG > Pintr, and Pintr > Pintr, suggesting there may be discrete (or non-linear manifolds of) synapse strength configurations that support bumps.

Prediction of Circuit Manipulation Effects
This framework allowed us to predict the effect of thermogenetic or optogenetic perturbation of neural populations. We computationally injected varying amounts of current into each neural subtype as defined by Wolff et al. (2015), i.e., distinguishing between E-PG and Pintr subtypes (Figure 5), and projected the ensuing circuit dynamics into the same space where we defined the dynamics modes ( Figure 4B). In general, the predicted effects matched the effects of changing the corresponding synapse class parameters. For example, injecting depolarizing current into the P-ENs had the same effect as increasing the strength of the excitatory E-PG > P-EN synapse class (or decreasing the strength of the inhibitory Pintr > P-EN synapse class). Injecting even relatively large (±5 nA) currents into the gall-tip-projecting subset of E-PGs or the P 6-8 -P 9 subset of Pintrs (Wolff et al., 2015) had little effect, presumably because these neural subtypes are less numerous in our model, represented by only two neurons each.

DISCUSSION
Ring attractor networks are an appealing explanation for the storage and updating of continuous variables in the brain (Skaggs et al., 1995;Zhang, 1996;Taube, 2007;Knierim and Zhang, 2012) and may play a role in visual attention (de Bivort and van Swinderen, 2016). We have shown dynamics consistent with a ring attractor arise in a network of generic spiking neurons with connectivity inferred from light-resolution microscopy and few other assumptions. The neurons in this network represent classes of neurons that are morphologically identical down to the level of independent neuropil subcompartments (glomeruli/wedges) defined in recent efforts to catalog all neurons in the PB of the CX (Wolff et al., 2015).
The model produces a number of key behaviors that are predicted by ring attractor theory (Song and Wang, 2005; FIGURE 5 | Prediction of circuit dynamics after physiological manipulation of neuronal subtypes. Circuit dynamics projected into two dimensions using the same input stimulus, diagnostic intervals, subsampling and linear projection as Figure 4B. Labeled neuronal subtypes were "injected" with ectopic currents as might be brought about by thermogenetic or optogenetic manipulation. Black lines represent 10 different current sweep replicates, and the thick color-mapped line their average. E-PG tip refers to the subset of E-PGs that project to the Gall tip. Pintr (P 6-8 -9 ) refers to intrinsic PB neurons that project from glomerulus 6-8 to glomerulus 9 and Pintr (P ∆7 ) to neurons that tile the PB with boutons while projecting dendrites throughout the PB. Taube, 2007) and observed by Seelig and Jayaraman (2015) by Ca 2+ activity in the E-PG neurons. In particular, a broad bump of activity (about 90 • -120 • wide) tracks a simulated cue as it moves. When the cue is removed, the bump persists indefinitely, initially staying in the position of the now-missing cue. In the continued absence of external cues, the bump can ''spontaneously'' shift in position. These shifts are absent when there is no spike noise in the input neurons ( Figure 3B). Thus, stochastic fluctuations in circuit input can move the bump, but it is hard to attribute these shifts to specific input spikes, perhaps not surprising considering the interconnected topology of the circuit.
We found that the bump may slide or jump to novel cues, and chooses only a single cue if multiple competitors are provided simultaneously (occasionally spontaneously jumping between them). Furthermore, we found that even when there is a pause in the representation of the bump by action potentials, it will reappear in the same position, as suggested by brief periods of time in which Seelig and Jayaraman (2015) observe the Ca 2+ bump to disappear and reappear in the same position. These findings suggest the bump is stored in subthreshold voltages.
Interestingly, our model suggests that there are discrete positions in the network in which this bump of activity prefers to reside as it moves through the network. Whether this is true of the circuit in vivo is not yet known, but it has been reported that startled cockroaches turn and run at angles that are multimodally distributed (Domenici et al., 2008). The modes of these escape angles are separated by approximately 30 • , which is nearly matched by the 13 modes of bump position that we observed ( Figure 2E). Perhaps discretized bump position tendencies underlie this distribution of escape angles. The origin of these modes surely arises from the discrete segmentation of the PB into nine glomeruli and the EB into 8 tiles/16 wedges. The modes appear to be evenly spaced every 22.5 • between ±112.5 • . Outside that range, two modes appear at approximately ±152 • . The bump dwells in these modes 1-2 orders of magnitude less often than the other modes. Perhaps their spatial irregularity (compared to the other modes) is related to their diminished stability. Furthermore, no mode is present at ±180 • , probably due to the broken radial symmetry of the circuit in this position ( Figure 1C), which is associated with the PB being split into linear segments while the EB has circular topology. The distribution of bump speeds ( Figure 2F) is also consistent with theoretical work predicting distinct modes of bump motion: sliding between adjacent positions and jumping between non-adjacent positions (Zhang, 1996). In general, sliding motion is induced when a stimulus moves slowly (or stochastic fluctuations increase near to the bump's current position). Conversely, the bump jumps when the stimulus moves rapidly (or stochastic fluctuations increase far from the bump's current position).
Additionally, angular position vectors can be coded not only by which neurons are active, but also by which pairs of neurons have synchronous activity (Ratté et al., 2013). In our circuit, we found that neurons tended to fire synchronously (Figure 3E), indicating that perhaps the PB could conceivably participate in a synchrony-based code. However, the sets of neurons observed to fire synchronously, particularly the P-EGs and P-ENs ( Figure 3A) are predicted by the number of monosynaptic connections each of these neuron classes is from the extrinsic inputs driving the circuit ( Figure 1D). Thus, the synchrony in the system may be determined by its topology and not flexible enough to usefully code angular position.
By adding random noise to synapse classes, we determined that a large range of synaptic strength parameters can result in apparently proper ring attractor dynamics (Figure 4). Moderate levels of noise in the synaptic strengths within the circuit often still produced dynamics consistent with experimental observations. We were able to characterize potential failure modes of the network, which include absence of action potentials in classes of neurons, an inability to sustain the bump, low responsiveness, low or no competitor selectivity, and/or network seizures.
By systematically varying the values of synapse strength parameters one at a time, we explored the space of failure modes to evaluate the robustness of our model. Our model predicts that perturbing certain synaptic or neuronal classes could have larger impacts on this network than others. In general, neuronal classes with fewer neurons could be perturbed more dramatically before causing a breakdown in bump dynamics. At the same time, all the synaptic or neuronal classes could be dramatically perturbed (or even reversed in sign) and still produce a proper bump, provided appropriate compensatory changes in other classes were made. Going forward, our model may be able to provide a quantitative framework for understanding variability in individual differences in navigation, such as locomotor handedness, a behavior whose variability across individuals is controlled by neurons in the PB (Ayroles et al., 2015;Buchanan et al., 2015). Specifically, we can add noise to individual synapses, reflecting variation from developmental stochasticity and post-developmental plasticity, e.g., by analogy to our synapse class parameter dithering: It is important to consider which assumptions made in this model might not be realistic. The information flow of each neuron class is inferred from the overlap of ''dendritic'' and ''axonal'' cellular compartments determined by lightmicroscopy. Despite being unipolar, neurons in Drosophila generally have polarized information flow (Rolls, 2011), however, common axo-axonal, dendro-dendritic and perhaps even dendro-axonal synapses (Schneider-Mizell et al., 2016) paint a more complex picture. Electrical coupling, which can lead to synchronized neuronal firing, is also common in insect neurons (Pereda, 2014), but we have not included any in our model. Furthermore, we have assumed that every neuron has the same integration and firing dynamics despite the fact that the dynamics can vary significantly based on specific ion channel expression levels (Marder, 2011;Berger and Crook, 2015). We also make the assumption that if an axon and dendrite overlap in a compartment then they are connected, but this is not necessarily the case. Neurons that are adjacent with the resolution provided by light microscopy may not come into physical contact (Feinberg et al., 2008). Moreover, axons and dendrites which are in contact do not necessarily form functional synapses (Kasthuri et al., 2015). Due to these caveats, it is remarkable that our model recapitulates so many of the experimental observations of Ca 2+ of E-PG neurons.
The core computation of this circuit may be robust to many categories of biological detail, emerging instead from high-level connectivity of the sort that can be inferred from light microscopy. Indeed, leaky-integrator implementations of the model simulating graded neurons without action potentials produced passable bumps (Figures 3F,G). This suggests that bump dynamics are not only robust to synapse strength parameter variation, but also variation in neuron physiological parameters. However, the leaky-integrator implementation did not exhibit selectivity between competitor stimuli as sharp as the leaky-integrate-and-fire implementation. This suggests that spiking neurons are required to recapitulate experimental Ca 2+ recordings, while not being required for bump dynamics in general.
Despite the conspicuous ring shape of the EB and its large number of inhibitory neurons with horizontal morphologies spanning all azimuthal positions (Martín-Peña et al., 2014;Kottler et al., 2017), neither of these qualities is necessary to bring about ring attractor-like dynamics in our model. Instead, our model generates global inhibition using intrinsic PB neurons (the PBintrs; Figure 1). The EB has been shown to receive spatiotopic information about visual features (Seelig and Jayaraman, 2013) and is involved in visual place learning (Ofstad et al., 2011). These observations suggest that the EB encodes spatial information about landmarks in the environment which could be used to correct accumulated error in the position of a bump. While inhibitory circuitry within the EB is not required for ring attractor dynamics in the PB-EB circuit, we have no evidence that the inhibitory circuitry in the EB does not participate in a ring attractor. It is possible that both the Pintrs in the PB and the ring neurons (Martín-Peña et al., 2014) of the EB implement long-range inhibition for the production of two distinct ring attractors, which could potentially interact to perform more sophisticated computations.
The egocentric heading correlate present in the PB-EB circuit is likely transmitted to other regions of the CX, particularly the FB. This neuropil could be a site of integration of navigational with internal state and sensory information for adaptive decision making. In addition to the PB-EB circuit neurons described here, the PB contains many columnar neurons projecting into the FB that have postsynaptic arbors in individual PB glomeruli and presynaptic boutons in different layers and columns of the FB (Wolff et al., 2015). Thus, it is likely that the FB inherits a bump or vertical band of activity from the PB. The FB is hypothesized to gate the selection of different behaviors in a state-dependent fashion (Weir and Dickinson, 2015) and activation of a single side of the FB induces ipsilateral turning (Guo and Ritzmann, 2013). Horizontal dopaminergic neurons in the FB have been shown to mediate sleep and arousal (Pimentel et al., 2016). The FB receives direct horizontal input from the visual system via the optic glomeruli (Mu et al., 2012) and also from many known modulatory neuropeptidergic neurons (Kahsai and Winther, 2011;Kahsai et al., 2012). The columnar projection neurons coming from the FB likely interact with these horizontal modulatory neurons. Therefore, it is appealing to hypothesize that the FB contains its own bump, downstream of the PB-EB bump, that it uses to integrate navigational information with neuromodulatory signals encoding internal states and sensory inputs.

AUTHOR CONTRIBUTIONS
BLB conceived of the project. KSK and BLB implemented the model, collected the data, analyzed it and wrote the manuscript.