Experimental setupFly preparation for imaging

We expressed the genetically encoded calcium indicator GCaMP7f (ref. 56) in EPG neurons by crossing GCaMP7f flies (w1118;;PBac[20XUAS-IVS-Syn21-op1-GCaMP7f-p10] in VK00005) to the EPG GAL4 driver line SS00096 (ref. 57). Flies (females, age 5–9 days, n = 10) were prepared for imaging as previously described8,58. Briefly, flies were anesthetized at 4 °C, their proboscis immobilized with wax to reduce brain movements, and their head/thorax fixed to a holder with a recording chamber using ultraviolet glue. To gain optical access to the brain, we removed a section of cuticle between the ocelli and antennae, along with the underlying fat and air sacs. Throughout the experiment, the head was submerged in saline containing NaCl (103 mM), KCl (3 mM), TES (5 mM), trehalose (8 mM), glucose (10 mM), NaHCO3 (26 mM), NaH2PO4 (1 mM), CaCl2 (2.5 mM) and MgCl2 (4 mM), with a pH of 7.3 and an osmolarity of 280 mOsm.

Two-photon calcium imaging

Calcium imaging was performed with a custom-built two-photon microscope controlled with ScanImage (version 2022, Vidrio Technologies)59. Excitation of GCaMP7f was generated with an infrared (920 nm), femtosecond-pulsed (pulse width ~110 fs) laser (Chameleon Ultra II, Coherent) with 15 mW of power, as measured after the objective (×60 Olympus LUMPlanFL/IR, 0.9 numerical aperture). Fast Z-stacks (eight planes with 6-μm spacing and three fly-back frames) were collected at 10 Hz by raster scanning (128 × 128 pixels, ~75 × 75 μm2) using an 8-kHz resonant-galvo system and piezo-controlled Z positioning. Focal planes were selected to cover the full extent of EPG processes in the EB. Emitted light was directed (primary dichroic: 735, secondary dichroic: 594), filtered (filter A: 680 SP, filter B: 514/44) and detected with a GaAsP photomultiplier tube (H10770PB-40, Hamamatsu).

Spherical treadmill system

Following dissection, flies were positioned on an air-supported polyurethane foam ball (8-mm diameter, 47 mg) under the two-photon microscope and allowed to walk. Rotations of the ball were tracked at 500 Hz, as described previously58. Behavioral data and imaging timestamps were recorded using WaveSurfer (version 0.947, http://wavesurfer.janelia.org/). For each fly, we collected five 20-min trials during which flies walked or stood in darkness.

Data analysis

All data analysis was performed in MATLAB (version 2022a, MathWorks). Some analyses relied on functions from the Circular Statistics Toolbox (version 2012a)60. No statistical methods were used to predetermine sample sizes, but our sample sizes are similar to those reported in previous publications8,30,61. Flies were selected at random from their vials; however, as all data were collected from a single experimental condition (flies walking in darkness), no other randomization was performed. Data collection and analysis were not performed blind to the conditions of the experiments. We excluded any data collected beyond 100 min for consistency and to exclude a small number of flies whose behavior and/or imaging degraded in quality, a known limitation of fly-on-a-ball calcium imaging experiments.

Extracting bump orientation and strength

Each Z-stack was reduced to a single frame using a maximum-intensity projection technique. An ellipse was manually drawn around the perimeter of the EB and automatically segmented into 32 equal-area, wedge-shaped ROIs. The number of ROIs was chosen to be twice the number of anatomically defined EB wedges62. Activity within each ROI was averaged for each frame, producing 32 ROI time series. For each ROI time series, baseline fluorescence (F0) was defined as the average of the lowest 10% of samples. ΔF/F was computed as 100 × (F − F0)/(F0), where F is the instantaneous fluorescence from the raw ROI time series. These ROI time series were then smoothed with a third-order Savitzky–Golay filter over 11 frames as in previous studies8,30. We used the PVA as a measure of bump strength and orientation. PVA was computed by taking the circular mean of vectors whose angles were the ROI’s wedge positions and whose length was equal to the ROI’s ΔF/F. The magnitude of this mean resultant vector length was normalized to have a maximum possible length of 1.

Characterizing bump drift

To determine bump drift (Fig. 1h,i), we first identified periods when flies were standing still (defined as zero rotational and translational velocity), disregarding periods shorter than 300 ms. Drift was computed as the circular distance between bump orientations (PVA phase) at the beginning and end of these periods of standing. To determine whether the EPG bump drifted from its initial position to preferred discrete locations within the EB when the fly stood still, we compared the distributions of initial and final bump positions across 64 nonoverlapping bins from −π to π around the structure (Extended Data Fig. 1a,b). We used Watson’s U2 test63,64, a nonparametric two-sample test, for this comparison, implemented using MATLAB code from P. Mégevand (watsons_u2, https://github.com/pierremegevand/watsons_u2, 2017). We used 500 permutations to compute P values for this test; these P values, together with the test statistic U2, are reported in the caption of Extended Data Fig. 1b. Finally, we computed the distribution of drifts for periods between 300 ms and 2 s across 64 nonoverlapping binned initial positions from −π to π around the EB, and fit each fly’s drift distribution with sinusoidal functions of the form A × sin(ω × ψ + θ) + C, where ω ∈  is the frequency of the sinusoid, ψ is the initial bump position during the standing period, and A, θ, C are learned parameters for the amplitude, phase, and DC offset, respectively (Extended Data Fig. 1c,d). Frequencies of 8 and 16 Hz were chosen to match the number of computational units in the fly’s compass network, which, in a discrete network, would cause the bump to drift toward 8 (or 16) distinct bump positions (schematized in Fig. 1h, top). For each fly, we computed the R2 value between the drift, measured as a function of HD, and the sinusoidal fits (Extended Data Fig. 1c); these R2 values are reported in Extended Data Fig. 1d.

Characterizing bump velocity

To determine whether the EPG bump shows signs of nonlinear integration (Fig. 1j, top), we measured whether the bump moved faster or slower than expected as a function of bump position for both left and right turns (Fig. 1j, middle and bottom). We began by performing a linear regression (ordinary least squares) between the fly’s instantaneous angular velocity and the bump’s angular velocity (both sampled at 10 Hz) to account for fly-to-fly variability in the gain of angular integration, as observed in previous studies8,30,61. Linear fits were separately performed for left and right turns, and the residuals were taken as a measure of whether the bump was moving faster (or slower) than expected after accounting for each fly’s naive gain. Next, we binned data by bump position (64 nonoverlapping bins from −π to π) and computed the average residual bump velocity for each bin, producing the curves shown in the middle and bottom panels of Fig. 1j. Lastly, we fit each fly’s curve with sinusoidal functions of the form A × sin(ω × ψ + θ) + C, where ω ∈  is the frequency of the sinusoid, ψ is the bump position, and A, θ, C are learned parameters for the amplitude, phase, and DC offset, respectively (Extended Data Fig. 2). Frequencies of 8 and 16 Hz were chosen to match the number of computational units in the fly’s compass network, which, in a discrete network, would cause the bump to move faster or slower than expected at 8 (or 16) distinct bump positions (schematized in Fig. 1j, top). For each fly, we computed the R2 value between the residual bump velocity, measured as a function of HD, and the sinusoidal fits (Extended Data Fig. 2a); these R2 values are reported in Extended Data Fig. 2b.

We note that our fly-on-a-ball calcium imaging setup comes with potential challenges for evaluating the presence or extent of nonlinear integration, including slow GCaMP dynamics, altered proprioceptive feedback that the fly may experience while walking on a ball heavier than itself, head fixation that may prevent the fly from altering its head–body angle during turns, potential neural propagation delays involved in relaying and integrating the angular velocity signal, and measurement noise inherent to calcium imaging that could corrupt bump velocity estimation.

Model overviewNetwork equations

We consider an effective single-ring network of N neurons (or, equivalently, of N computational units; see ‘Network equations’ in the Supplementary Note). Neurons are ordered according to their preferred heading θj, which we take to be evenly spaced by Δθ = 2π/N rad. Neurons are recurrently connected according to their preferred headings through a symmetric weight matrix \(_^}}=_}+_}\cos (_-_)\), where JE and JI parametrize the strength of local excitation and uniform inhibition, respectively (note that JE and JI actually correspond to tuned and untuned components of the connectivity; for ease of language, we use local excitation and broad inhibition here and throughout). Neurons receive velocity input through an asymmetric, velocity-modulated weight matrix \(_}}_^}}=_}}\sin (_-_)\); in the main text, we took vin > 0. Each neuron j receives a constant feedforward input cff and a net input \(1/N\,_(_^}}+_}}_^}})_\) from all other neurons in the network, where the firing rate rk = ϕ(hk) is a nonlinear function of the total input activity hk. For all analyses shown in the main text, we took the nonlinear transfer function ϕ(⋅) to be rectified linear (that is, ϕ(⋅) = [⋅]+, but see also Extended Data Fig. 5 and ‘Robustness to changes in the transfer function and recurrent weights’ in the Methods). The dynamics of each neuron are given by the system of single-neuron equations in equation (1); we chose τ = 0.1 s and cff = 1.

By applying a discrete Fourier transform to the single-neuron equations, we can express this system of equations in terms of its Fourier modes. After initial transients, only the DC and first-order modes remain, and the resulting dynamical system reduces to a set of three equations that govern the dynamics of the orientation ψ, amplitude a relative to the average input activity, and width w of the bump (‘Order equations’ in the Supplementary Note); we will refer to these as the system of bump equations.

Stable parameter regime

The system of bump equations will generate a stable bump of activity for certain combinations of JE and JI (‘Fixed point analysis’ in the Supplementary Note and Extended Data Fig. 3a). For all analyses shown in the main text, we first selected a desired value of JE > 2, and then selected a value of JI such that it produced a bump of activity whose full amplitude A = H0 + a (where H0 is the average input activity) was at least approximately A = 0.2. To do so, we first uniformly sampled bump orientations ψ ∈ [0, 2π) and widths w ∈ [2π/N, 2(N − 1)π/N), and we used these to calculate the contour JEfeven(w, ψ) = 1 using MATLAB’s ‘contourc.m’, where feven(w, ψ) is given by equation (S19) in the Supplementary Note (see also equation (S30) in the Supplementary Note and Extended Data Fig. 3c). This gave us values \((w,\psi )\in __}}=\_}_}}(w,\psi )=1\}\) that satisfy the contour equation. We then used these values of w and ψ to determine an upper bound on JI given by

$$_}^}}=\mathop\limits___}}}\frac_(w,\psi )},$$

(8)

where f0(w, ψ) is given by equation (S18) (see also equation (S32)) in the Supplementary Note. We then used these same values of w and ψ to determine a value for JI, given by

$$_}=\mathop\limits___}}}\frac_}/A-1)\cos (w/2)-_}/A}_(w,\psi )},$$

(9)

and verified that \(_} < _}^}}\). Plugging A = 0.2 into equation (9) resulted in a bump of activity whose minimum full amplitude was approximately A = 0.2.

Model analyticsStationary solutions

To determine the configurations to which the system evolves in the absence of velocity input, we characterized the stationary solutions of the system of bump equations (‘Fixed point analysis’ in the Supplementary Note). This allowed us to determine relationships between the bump orientation, relative amplitude, and width that would persistently maintain a stable bump of activity (Extended Data Fig. 3b,c). For a network of N neurons that receive no velocity input, most parameter settings will yield two sets of N fixed points each—one set will be stable, and the other will be unstable. For a given value of JE, one set will be aligned with the preferred headings , and the other set will be aligned precisely between the preferred headings; the second and fourth columns of Fig. 2e highlight examples for which the unstable (second column) and stable (fourth column) sets of fixed points are aligned with the preferred headings. The value of JE and the parity of N (whether the network consists of an even or odd number of neurons) together specify which of these two configurations the network will adopt. When N is even and \(_} < _},N-2}^\) (denoting bumps supported by N − 1 and N − 2 neurons), the set of fixed points aligned with the preferred headings will be unstable. When N is odd, the reverse will be true: for \(_} < _},N-2}^\), the set of fixed points aligned with the preferred headings will be stable. For a given network size N, as JE passes through an optimal value \(_}^\), this stability switches (Extended Data Fig. 3d,g). At each of these fixed points, the widths of the stable and unstable bump configurations are determined solely by JE, whereas their relative amplitudes depend on both JE and JI.

Energy landscape

We derived an energy landscape E(a, w, ψ; JE, JI) for the system of bump equations in the absence of velocity input40,41 (‘Energy landscape’ in the Supplementary Note). This function describes the stable configurations to which the system will evolve in the absence of input.

To minimize the curvature of the energy landscape, we first determined the 3 × 3 Hessian matrix of the second derivatives of the energy E with respect to a, w, and ψ. When evaluated at the orientations ψs of the stable fixed points (see the previous subsection), we found that the Hessian reduced to a block diagonal matrix, with a single eigenvector along ψ whose eigenvalue is given by

$$\frac^E}^}\propto 1-\frac_}}\sum __}}}^(_-^}),$$

(10)

where Kact denotes the set of indices of the neurons that actively maintain the bump. This eigenvalue quantifies the degree of local curvature of the energy as a function of bump orientation ψ. For a system of size N, there are N − 3 values of local excitation JE for which this eigenvalue goes to zero, and thus for which the energy landscape is locally flat as a function of ψ. These correspond to bump configurations for which the bump is maintained by Nact ∈ [2, N − 2] active neurons:

$$\frac_},_}}}^}=\frac+\frac\left(\tilde+\frac/N)}\right);\,\tilde=_}}-\frac.$$

(11)

We found that these values of local excitation, which are shown in Fig. 2d, also ensure that the energy landscape is flat for all bump orientations (as shown in Fig. 2e; also see Extended Data Fig. 4).

Leading eigenvalues of active submatrices

In the absence of velocity input, the bump dynamics are governed by the leading eigenvalue λ of a submatrix of the connectivity (−I + Wsym/N)/τ; this eigenvalue determines the rate at which the bump will drift in the absence of input. When the local excitation JE is optimally tuned (that is, \(_}=_},_}}}^\)), the bump of activity will be maintained by a fixed number of active neurons Nact ∈ [2, …, N − 2]. For each distinct value of Nact, there is thus a distinct Nact × Nact submatrix of the connectivity whose single leading eigenvalue determines the drift dynamics. Away from these optimal values of local excitation, the bump of activity will be maintained by either n or n + 1 active neurons (see equation (S50) in the Supplementary Note). The drift dynamics are then governed by the leading eigenvalues of the corresponding n × n and (n + 1) × (n + 1) active submatrices.

To determine these dynamics, we analytically determined the rates of bump drift in the stable and unstable regimes, which are given in equation (2) (see ‘Performance of non-optimal solutions: Dynamics in the absence of input velocity’, and, in particular, equations (S54) and (S56) in the Supplementary Note). We then compared these analytically derived drift rates to the leading eigenvalues that we computed numerically by directly diagonalizing active submatrices of the connectivity (using the MATLAB function ‘eig.m’); this comparison is shown in Extended Data Fig. 6.

Widths of stable and unstable regimes

In the absence of input, the widths of the stable and unstable regimes can be determined analytically by finding the orientation at which the bump transitions from unstable to stable dynamics as it drifts away from an unstable fixed point. This reduces to matching two exponential equations that govern the dynamics of the bump orientation in the two regimes (with drift rates λu and λs, respectively), and that must tend toward the orientations of the unstable and stable fixed points as t → −∞ and t → +∞, respectively. The resulting widths of each regime are given by equation (3) and shown in Fig. 3d and Extended Data Fig. 8b, and they are centered on the orientations of the stable and unstable fixed points in the absence of input. Given a stable fixed point at ψ = ψs and an unstable fixed point at ψ = ψu = ψs + π/N, the resulting equation for the bump can then be written as (see equations (S61) and (S62) in the Supplementary Note):

$$\psi (t)=\left\^}+(_-^})\exp (_}t)\,\,;\quad 0 < t < _\quad \,\text\,\quad \\ ^}+\frac_}}\exp (-| _}| (t-_))\,;\quad t > _\quad \,\text\,\quad \end\right.,$$

(12)

where ψs + Δθs/2 < ψ0 < ψu is the initial orientation of the bump, and tΔn = (1/λu) log(Δθu/(2(ψu − ψ0))) is the time when the bump orientation crosses from the unstable regime into the stable regime. See ‘Performance of non-optimal solutions: Dynamics in the absence of input velocity’ in the Supplementary Note for more details.

Drift in the absence of input

To measure the net bump drift, we analytically computed the time τd that it takes for the bump to drift from within εu of an unstable fixed point to within εs of a stable one. We chose εu = Δθu/2e and εs = Δθs/2e, such that the bump covered an angular distance of Δψd = (1 − 1/e)Δθ/2 in the time τd. We then measured the net drift speed as Δψd/τd (see equations (S68)–(S71) in the Supplementary Note).

Small velocity approximation

In the presence of velocity input, the bump dynamics will be governed by the leading eigenvalue λ of a submatrix of the full connectivity (−I + (Wsym + vinWasym)/N)/τ. The asymmetric component of this connectivity is modulated by the input velocity vin, and introduces a velocity-dependent correction to the eigenvalue λ0 of the symmetric connectivity (−I + Wsym/N)/τ (Extended Data Fig. 7):

$$\lambda \approxeq _+f(\,_})_}}^+}(_}}^).$$

(13)

For sufficiently small input velocities, we can approximate the leading eigenvalues λu and λs, and thus the corresponding widths of the unstable and stable regimes, as being equal to their values in the absence of velocity input (see ‘Leading eigenvalues of active submatrices’ and ‘Widths of stable and unstable regimes’ in the Methods). All analytic results shown in Fig. 3i–l were generated under this assumption. This approximation breaks down as the input velocity increases, and it breaks down more quickly for smaller values of local excitation (as shown in Fig. 3l; see also Extended Data Fig. 7a).

Locations of fixed points in a velocity-driven regime

Although we can approximate the rates and width of the stable and unstable regimes as remaining unchanged for a sufficiently small velocity input, we cannot make the same approximation for the orientations of stable and unstable fixed points. Therefore, we will treat the stable and unstable fixed-point orientations as functions of vin: ψs = ψs(vin), ψu = ψu(vin), respectively. The orientation of the stable and unstable fixed points found in the absence of velocity input will then be given by ψs(0) and ψu(0), respectively. To determine how the orientations of these fixed points shift with velocity, we repeated the analyses described in ‘Widths of stable and unstable regimes’ in the Methods, but with a different set of initial conditions (see ‘Performance of non-optimal solutions: Dynamics in the presence of small input velocity’ in the Supplementary Note for details). Given a bump that begins at a stable fixed point ψ = ψs(0) in the absence of input, and given an initial velocity vin, the bump will be driven to a new stable fixed point at an orientation ψs(vin) = ψs(0) + vin/|λs| as t → ∞. In the limit that t → −∞, the bump will be driven to (and hence, in forward time, away from) an unstable fixed point at an orientation ψu(vin) = ψu(0) − vin/λu. Over an interval ψ ∈ [ψs(0) − Δθs/2, ψu(0) + Δθu/2], the resulting equation for the bump can be written as (see equations (S78) and (S79) in the Supplementary Note):

$$\begin\psi (t)=\left\^}(0)+\frac_}}}_}| }\left(1-\exp (-| _}| t )\right);\qquad \qquad0 < t < _}\;\,\text\\ ^}(0)-\frac_}}}_}}+\left(\frac_}}}_}}-\frac(\Delta \theta -\Delta _})\right)\\\qquad\;\;\exp (_}(t-_}));\;\qquad\qquad\qquad \quad\;\;> _}\qquad\;\text \end\right.,\end$$

(14)

where tc = (1/|λs|) log(1/(1 − Δθs|λs|/2vin)) is the time when the bump orientation crosses from the stable regime into the unstable regime.

At the threshold velocity given in equation (5), the two fixed points will meet at the boundary between regimes; this is the minimum velocity needed for the bump to move continuously. Below this velocity, the bump will be driven away from the unstable fixed point in the unstable regime, and toward a stable fixed point in the stable regime. Above this velocity, the stable and unstable fixed points will still drive the bump dynamics, but their orientations will move outside of their respective regimes. The minimum and maximum bump velocities, \(_}\) and \(_}\) (given by equation (6)), can be computed analytically from equation (14) by evaluating the time derivative of ψ(t) at the boundary from the stable to the unstable regime, and vice versa. We used these minimum and maximum velocities to define the linearity of integration as \(_}/_}\). See ‘Performance of non-optimal solutions: Dynamics in the presence of small input velocity’ in the Supplementary Note for details.

Simplified energy landscape

Having described each linear subsystem in terms of (1) the orientations of the fixed points, (2) the rate at which the bump drifts toward or away from these fixed points, and (3) the angular regime governed by each fixed point, we used these three properties to construct a simplified landscape that describes the energy of different bump orientations. Given a linear system, an energy function can be chosen to be quadratic65; we thus choose Eu,s(ψ) = αu,sψ2, where αs > 0 for the stable subsystem, and αu < 0 for the unstable subsystem. To select the appropriate values of αu,s, we require that the energy function has extrema at the orientations of the stable and unstable fixed points ψs(vin) and ψu(vin), and that the energy transitions smoothly between the stable and unstable regimes; this yields

$$E(\psi )=\left\-_}^}(_}}))}^+C(_}})\quad \,\text\,\quad \\ -_}^}(_}}))}^\quad \quad \quad \quad \quad \,\text\,\quad \end\right.,$$

(15)

where \(C(_}})=\Delta \theta \left(_}\Delta _}/4-_}}+_}}^/(_}\Delta _})\right)\), and where αu = −λu < 0 and αs = −λs > 0 as required. When moving around the ring, each successive pair of stable and unstable regimes will be governed by an energy landscape of this form but with a vertical shift, such that E(ψ ± nΔθ) = E(ψ) ∓ 2nvinΔθ. See ‘Simplified energy’ in the Supplementary Note for more details.

Tolerance in tuning

To determine how precisely the local excitation must be tuned to achieve a criterion level of performance, we first computed the derivative of each performance measure as a function of local excitation, evaluated at an optimal value; we denote this \(_(\,_}^)\) (see equations (S96)–(S99) in the Supplementary Note). This slope gives us a local linear estimate of how quickly the performance degrades away from an optimal value of local excitation. Because each performance measure can be expressed as a function of the net drift speed |λd|, computing this slope reduced to computing \(\partial | _}| /\partial _}__}^}\). Given a criterion for the system to be within \(_^}}\) of optimal performance for a performance measure P, the tolerance about a given optimal value \(_}^\) can then be computed as \(\left.}_(\,_}^)\ge 2_^}}\right/| _(\,_}^)|\) (where ≥ indicates that this is a lower bound on the tolerance, as the linear slope will overestimate the rate of degradation of performance; see equation (S113) in the Supplementary Note).

To determine the volume of parameter space that can meet this desired performance, we summed the tolerance across all optimal values of local excitation for a given network size N (see equation (S120) in the Supplementary Note). We then approximated this sum by its largest value, which reduces to

$$V(N)\ge _\frac^}.$$

(16)

See ‘Degradation of performance as a function of local excitation’ in the Supplementary Note for more details.

Model simulationsOverview

All simulations that we performed used MATLAB’s ODE solver ‘ode45.m’ with an integration timestep of Δt = 0.01 s. We first initialized the network to generate a bump of activity at a given orientation ψ. Using this as the initial condition for the network, we then simulated the single-neuron dynamics in equation (1), and we performed a discrete Fourier transform using MATLAB’s ‘fft.m’ function to extract the bump dynamics as a function of the single-neuron dynamics (see equation (S16) in the Supplementary Note). When simulating angular velocity integration, we first determined the velocity scaling that would generate a comparable rate of bump movement for a given (constant) velocity input (see ‘Velocity-driven dynamics’ in the Methods). We then simulated the network dynamics in response to this scaled input.

Parameter choices

All results shown in Figs. 2 and 3 were generated using networks of size N = 6. When illustrating network properties for different values of local excitation, we used the following values of JE (evenly spaced in 1/JE): JE = [12, 6, 4, 3, 2.4] (Fig. 2e–h); JE = [3.89, 3.6, 3.3, 3, 2.77, 2.57, 2.44] (Fig. 3f,j); JE = [3.6, 3, 2.57, 2.44] (Fig. 3g,k); 17 evenly spaced values of 1/JE between 1/3.86 and 1/2.45 (Fig. 3h,l). When simulating network dynamics in the presence of velocity input, we used the following values of velocity input vin: ten evenly spaced velocity values between 0.2 and 2.0 rad s−1 (Fig. 2f); ten evenly spaced values between 0.1 and 1.0 rad s−1 (Fig. 3k); five evenly spaced values between 0.8 and 1.6 rad s−1 (Fig. 3l). In all cases, we scaled the velocity input as described below (see ‘Velocity-driven dynamics’ in the Methods).

Drift in the absence of input

For simulations of bump drift, we simulated the network with the velocity input set to zero. To illustrate drift trajectories for different values of JE (as shown in the bottom row of Fig. 2f and in Fig. 3g), we initialized the bump at six evenly spaced orientations between (and including) 0 and π/N, and we simulated the evolution of the bump for 3 s. We repeated this for repeating angular units between 0 and 2π.

Measuring net drift speed

To measure the net drift speed (as described in ‘Drift in the absence of input’ in ‘Model analytics’ in the Methods), we initialized the bump at an orientation ψu − εu (where ψu is the orientation of an unstable fixed point; for the values of JE used in Fig. 3, ψu = π/N; see ‘Parameter choices’ in the Methods). We then simulated the network dynamics until the bump reached an orientation εs. We set εu = Δθu/2e and εs = Δθs/2e, where Δθu,s were computed as described in ‘Widths of stable and unstable regimes’ in the Methods. We used the time it took for the bump to reach this orientation as the measure of the net drift timescale τd, and we used Δψd/τd as a measure of net drift speed, where Δψd = (1 − 1/e)Δθ/2 is the angular distance traveled by the bump in the time τd. Fig. 3h compares the net drift speed from simulations to that obtained analytically for different values of JE.

Velocity-driven dynamics

For simulations of angular velocity integration, we injected a constant velocity input throughout the simulation. To permit a comparison to analytic predictions, we scaled the input velocity such that the rate of movement of the bump matched the input velocity at an input of vin = 50 rad s−1. To this end, we determined the best-fitting linear trajectory that minimized the absolute deviation from the bump trajectory over a time window of t = 6 s, and we used the slope of this linear trajectory to scale all other input velocities injected into the network. We performed this scaling separately for each set of network parameters (that is, for each choice of (JE, JI)). All velocity values described in simulations were scaled in this way.

Measuring threshold velocity

To measure the threshold velocity required to move the bump continuously (as shown in Fig. 3l), we first analytically computed the threshold velocity as described in ‘Locations of fixed points in a velocity-driven regime’ in the Methods. We then chose 50 evenly spaced input velocity values between (and including) vthresh − 0.05 rad s−1 and vthresh + 0.05 rad s−1. We initialized the bump at the orientation of a stable fixed point (here, at ψs = 0), and we then simulated the network dynamics in response to each velocity individually. We determined the minimum of these velocities that would move the bump beyond an orientation of π/N within a time interval of 10 s. Fig. 3l compares this simulated value to the value obtained analytically.

Measuring the linearity of integration

To measure the linearity of integration from simulations, we simulated the bump trajectory for different constant input velocities (as described above in ‘Overview’). For each input velocity, we determined the time tc when the bump orientation ψ crossed from the stable into the unstable regime or vice versa; these times were used to compute the minimum and maximum velocities, respectively (note that we used the analytically derived boundaries between regimes to determine these crossing times; see ‘Widths of stable and unstable regimes’ in the Methods). We then determined the bump velocity as ν = (ψ(tc + Δt) − ψ(tc − Δt))/2Δt, where Δt = 0.1 s is the integration timestep used in all simulations. Fig. 3l compares this simulated value to the value derived analytically (see ‘Locations of fixed points in a velocity-driven regime’ in the Methods).

Robustness to variations in parameter tuning

To summarize performance as a function of network size (shown in Fig. 4a), we analytically computed the net drift speed (as described in ‘Drift in the absence of input’ in ‘Model analytics’ in the Methods) as a function of local excitation in the range \(_}\in [\,_},N-2}^,_},2}^]\) (that is, between the minimum and maximum optimal values of local excitation, maintained by Nact = N − 2 and Nact = 2 active neurons, respectively). For each optimal value of local excitation, we numerically estimated the tolerance as the range of local excitation values about an optimum for which the net drift speed would be consistently below a fixed performance threshold (we used a threshold value of 0.001 rad s−1). We considered only those values of local excitation above the minimum optimal value or below the maximum optimal value to estimate this tolerance; thus, to estimate the tolerance about the minimum and maximum optimal values, we measured the tolerance in only one direction (\(_}\le _,2}^\) or \(_}\ge _},N-2}^\)), and we doubled this value to use as our estimate. We then compared these tolerance estimates to the analytic lower bound given in equation (7), as shown in Fig. 4b,c (also see equations (S113)–(S119) in the Supplementary Note). Finally, we summed these tolerance values (computed numerically or analytically) for each network size N to estimate the net volume of parameter space that meets this threshold level of performance, as shown in Extended Data Fig.